text:
#=start.jl - Starts the game and accepts 7 function call aruments:
=#
include("board.jl")
module start

using BM

using SQLite

function startShogi(database,gameTypein,cheatingin,timed,time_add,sente_time,gote_time)

#  database = ARGS[1] #/path/to/database/file {string}
  gameTypeDict = Dict("S" => "standard","M" => "minishogi")
  gameType = gameTypeDict[gameTypein] #Either "standard" or "minishogi" {string}
  cheatingDict = Dict("T" => "cheating", "F" => "legal")
  cheating = cheatingDict[cheatingin] #Either "cheating" or "legal" {string}=#
  #timed is "yes" or "no"
  #captures = Array(square(),0)
  seed = time() #current unix time
  db = SQLite.DB(database) #Opens the database gamefile
  SQLite.query(db,"CREATE TABLE moves (move_number,move_type,sourcex,sourcey,targetx,targety,option,i_am_cheating,targetx2,targety2,targetx3,targety3);")
  SQLite.query(db,"CREATE TABLE meta (key,value);")
  SQLite.query(db,"""INSERT INTO meta (key,value) VALUES ("type","$(gameType)");""")
  SQLite.query(db,"""INSERT INTO meta (key,value) VALUES ("legality","$(cheating)");""")
  SQLite.query(db,"""INSERT INTO meta (key,value) VALUES ("seed","$(seed)");""")
  SQLite.query(db,"""INSERT INTO meta (key,value) VALUES ("timed","$(timed)");""")
  SQLite.query(db,"""INSERT INTO meta (key,value) VALUES ("time_add","$(time_add)");""")
  SQLite.query(db,"""INSERT INTO meta (key,value) VALUES ("sente_time","$(sente_time)");""")
  SQLite.query(db,"""INSERT INTO meta (key,value) VALUES ("gote_time","$(gote_time)");""")

end
export startShogi
end

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia


text:
using Datetime

y = year(1)
m = months(1)
w = weeks(1)
d = days(1)
h = hour(1)
mi = minute(1)
s = second(1)

#Period arithmetic
Base.Test.@test typeof(y+m) ==  Month{ISOCalendar}
Base.Test.@test typeof(m+y) ==  Month{ISOCalendar}
Base.Test.@test typeof(y+w) ==  Week{ISOCalendar}
Base.Test.@test typeof(y+d) ==  Day{ISOCalendar}
Base.Test.@test typeof(y+h) ==  Hour{ISOCalendar}
Base.Test.@test typeof(y+mi) == Minute{ISOCalendar}
Base.Test.@test typeof(y+s) ==  Second{ISOCalendar}
Base.Test.@test_throws typeof(m+w)
Base.Test.@test_throws typeof(m+d)
Base.Test.@test_throws typeof(m+h)
Base.Test.@test_throws typeof(m+mi)
Base.Test.@test_throws typeof(m+s)
Base.Test.@test typeof(w+d) ==  Day{ISOCalendar}
Base.Test.@test typeof(w+h) ==  Hour{ISOCalendar}
Base.Test.@test typeof(w+mi) ==  Minute{ISOCalendar}
Base.Test.@test typeof(w+s) ==  Second{ISOCalendar}
Base.Test.@test typeof(d+h) ==  Hour{ISOCalendar}
Base.Test.@test typeof(d+mi) ==  Minute{ISOCalendar}
Base.Test.@test typeof(d+s) ==  Second{ISOCalendar}
Base.Test.@test typeof(h+mi) ==  Minute{ISOCalendar}
Base.Test.@test typeof(h+s) ==  Second{ISOCalendar}
Base.Test.@test typeof(mi+s) ==  Second{ISOCalendar}
Base.Test.@test y * m*months(2) == months(24)
Base.Test.@test y * w*weeks(2) == weeks(104)
Base.Test.@test y * d*days(2) == days(730)
Base.Test.@test y > m
Base.Test.@test d < w
Base.Test.@test typemax(Year) == 2147483647
Base.Test.@test typemax(Year) + y == -2147483648
Base.Test.@test typemin(Year) == -2147483648
q = convert(Year{Calendar},1)
Base.Test.@test typeof(q) == Year{Calendar}
y,q = promote(y,q) #check periods w/ diff calendars
Base.Test.@test typeof(q) == Year{ISOCalendar}
Base.Test.@test typeof(y) == Year{ISOCalendar}
Base.Test.@test q + y == years(2)
#Period-Real arithmetic
Base.Test.@test y + 1 == 2
Base.Test.@test 1 + y == 2
Base.Test.@test y + 1.0 == 2.0
Base.Test.@test y * 4 == 4
Base.Test.@test y * 4f0 == 4.0f0
Base.Test.@test y * 3//4 == year(1)
Base.Test.@test y / 2 == year(1)
Base.Test.@test 2 / y == 2
Base.Test.@test y / y == 1
Base.Test.@test y*10 % 5 == 0
Base.Test.@test 5 % y*10 == 0
Base.Test.@test !(y > 3)
Base.Test.@test 4 > y
t = [y,y,y,y,y]
Base.Test.@test t .+ year(2) == [year(3),year(3),year(3),year(3),year(3)]
#PeriodRange
t = year(100)
r = y:t
Base.Test.@test size(r) == (100,)
Base.Test.@test length(r) == 100
Base.Test.@test step(r) == year(1)
Base.Test.@test first(r) == y
Base.Test.@test last(r) == t
Base.Test.@test first(r + y) == year(2)
Base.Test.@test last(r + y) == year(101)
Base.Test.@test last(y:year(2):t) == year(99)
Base.Test.@test [year(0):year(25):t][4] == year(75)
Base.Test.@test_throws y:m:t

#traits
Base.Test.@test zero(typeof(y)) == year(0)
Base.Test.@test one(typeof(y)) == y
Base.Test.@test months(y) == months(12)
Base.Test.@test weeks(y) == weeks(52)
Base.Test.@test days(y) == days(365)
#wrapping
Base.Test.@test addwrap(1,11) == 12 #month 1 plus 11 months == month 12
Base.Test.@test addwrap(1,12) == 1  #month 1 plus 12 months == month 1
Base.Test.@test addwrap(1,24) == 1  #month 1 plus 24 months == month 1
Base.Test.@test subwrap(12,1) == 11 #month 12 minus 1 month == month 11
Base.Test.@test addwrap(2000,12,1) == 2001 #year 2000 month 12 plus 1 month == year 2001

#Date tests
dt = date(1,1,1)
Base.Test.@test convert(Date,1) == dt
Base.Test.@test convert(Date,int32(1)) == dt
Base.Test.@test convert(Date,float32(1)) == dt
dt2 = convert(Date{Calendar},float64(1))
Base.Test.@test typeof(dt2) == Date{Calendar}
a,b = promote(dt,dt2) #checks 2 diff calendars promotion
Base.Test.@test typeof(a) == Date{ISOCalendar}
Base.Test.@test typeof(b) == Date{ISOCalendar}
Base.Test.@test dt == dt2
Base.Test.@test dt - dt2 == days(0)
Base.Test.@test_throws dt + dt2
Base.Test.@test !(dt > dt2)
Base.Test.@test int64(hash(dt)) == 6614235796240398542
Base.Test.@test length(dt) == 1
Base.Test.@test convert(Date,0) == date(year(0),month(12),day(31))
Base.Test.@test convert(Date,-1) == date(0,12,30)
Base.Test.@test_throws dt2 + y #Diff calendars
dt3 = datetime(1,1,1)
Base.Test.@test convert(DateTime,86400000) == dt3
dt4 = convert(DateTime{Calendar,UTC},float64(86400000))
Base.Test.@test typeof(dt4) == DateTime{Calendar,UTC}
a,b = promote(dt3,dt4) #checks 2 diff calendars promotion
Base.Test.@test typeof(a) == DateTime{ISOCalendar,UTC}
Base.Test.@test typeof(b) == DateTime{ISOCalendar,UTC}
Base.Test.@test dt3 == dt4
Base.Test.@test dt3 - dt4 == seconds(0)
Base.Test.@test_throws dt3 + dt4
Base.Test.@test !(dt3 > dt4)
Base.Test.@test int64(hash(dt3)) == 2879920070571771443
Base.Test.@test length(dt3) == 1
a,b = promote(dt,dt3) #promote Date => DateTime
Base.Test.@test typeof(a) == DateTime{ISOCalendar,UTC}
Base.Test.@test typeof(b) == DateTime{ISOCalendar,UTC}
Base.Test.@test date(dt3) == dt
Base.Test.@test datetime(dt) == dt3
Base.Test.@test dt == dt3
Base.Test.@test !(dt > dt3)
Base.Test.@test dt - dt3 == seconds(0)
Base.Test.@test_throws dt + dt3

check = (-252522163911150,-1,0,1,252522163911149)
for (i,x) in enumerate(check)
	yd = Datetime.yeardays(x-1)
	md = Datetime.monthdays[1]
	d = 1
	td = Datetime.totaldays(x,1,1)
	Base.Test.@test d + md + yd - 306 == td
	Base.Test.@test td == (-92231826452318568,-730,-365,1,92231826452317474)[i]
	Base.Test.@test Datetime._year(td) == x
	Base.Test.@test Datetime._month(td) == 1
	Base.Test.@test Datetime._day(td) == 1
	Base.Test.@test Datetime._day2date(td) == (x,1,1)
	Base.Test.@test Datetime._week(td) == (1,1,1,1,53)[i]
end
Base.Test.@test isleap(-4)
Base.Test.@test isleap(0)
Base.Test.@test isleap(4)
Base.Test.@test !isleap(100)
Base.Test.@test isleap(400)
Base.Test.@test isleap(2000)
Base.Test.@test !isleap(2100)

dt = date(2013,7,1)
Base.Test.@test year(dt) == 2013
Base.Test.@test month(dt) == 7
Base.Test.@test day(dt) == 1
Base.Test.@test Datetime._days(dt) == 735050 == int64(dt)
Base.Test.@test string(dt) == "2013-07-01"
Base.Test.@test week(dt) == 27
Base.Test.@test isdate(dt)
Base.Test.@test calendar(dt) == ISOCalendar
Base.Test.@test typemax(Date) == date(252522163911149,12,31)
Base.Test.@test typemin(Date) == date(-252522163911150,1,1)
Base.Test.@test isleap(dt) == false
Base.Test.@test lastdayofmonth(dt) == 31
Base.Test.@test dayofweek(dt) == 1
Base.Test.@test dayofweekinmonth(dt) == 1
Base.Test.@test daysofweekinmonth(dt) == 5
Base.Test.@test firstdayofweek(dt) == dt
Base.Test.@test lastdayofweek(dt) == dt + days(6)
Base.Test.@test dayofyear(dt) == 182
Base.Test.@test typeof(today()) <: Date
dt = date(2012,2,29)
dt2 = date(2000,2,1)
Base.Test.@test dt > dt2
Base.Test.@test dt != dt2
Base.Test.@test year(-dt) == -2012
Base.Test.@test +dt == dt
Base.Test.@test_throws dt + dt2
Base.Test.@test_throws dt * dt2
Base.Test.@test_throws dt / dt2
Base.Test.@test_throws dt + hour(1)
Base.Test.@test dt - dt2 == 4411
Base.Test.@test dt > dt2
Base.Test.@test dt2 < dt
Base.Test.@test dt2 - year(3) == date(1997,2,1)
Base.Test.@test dt + year(1) == date(2013,2,28)
Base.Test.@test year(3) - dt2 == dt2 - year(3)
Base.Test.@test dt2 - months(3) == date(1999,11,1)
Base.Test.@test months(11) + dt == date(2013,1,29)
Base.Test.@test months(8) + dt == dt + months(8)
Base.Test.@test dt2 + days(4411) == date(2012,2,29)
Base.Test.@test dt2 + days(4412) == date(2012,3,1)
Base.Test.@test dt2 + weeks(52) == date(2001,1,30)
Base.Test.@test dt2 + weeks(104) == date(2002,1,29)
Base.Test.@test dt2 - year(4) + days(366) == date(1997,2,1)
Base.Test.@test dt2 - (year(4) - days(366)) == date(1997,2,2) #non-associative
dt = date(1999,12,27)
check = (52,52,52,52,52,52,52,1,1,1,1,1,1,1,2,2,2,2,2,2,2)
for i = 1:21
	Base.Test.@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = date(2000,12,25)
for i = 1:21
	Base.Test.@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = date(2030,12,23)
for i = 1:21
	Base.Test.@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = date(2004,12,20)
check = (52,52,52,52,52,52,52,53,53,53,53,53,53,53,1,1,1,1,1,1,1)
for i = 1:21
	Base.Test.@test week(dt) == check[i]
	dt = dt + day(1)
end
#Simple Date parsing
a = "2005-01-31"
b = "2005/31/01"
c = "01.31.2005"
d = "31 01 2005"
aa = "2005,1,31"
bb = "2005-31-1"
cc = "1-31-2005"
dd = "31-1-2005"
aaa = "05-01-31"
bbb = "05/31/01"
ccc = "01.31.05"
ddd = "31 01 05"
e = "20050131"
ee = "2005131"
for i in (a,b,c,d,aa,bb,cc,dd,ccc,ddd)
	Base.Test.@test date(i) == date(2005,1,31)
end
Base.Test.@test date(aaa) == date(2031,5,1)
Base.Test.@test date(bbb) == date(2001,5,31)
Base.Test.@test_throws date(e)
Base.Test.@test_throws date(ee)

#DateRange
dt1 = date(2000,1,1)
dt2 = date(2010,1,1)
y = year(1)
r = dt1:y:dt2
Base.Test.@test first(r) == date(2000,1,1)
Base.Test.@test last(r) == date(2010,1,1)
Base.Test.@test typeof(r.step) == Year{ISOCalendar}
Base.Test.@test length(r) == 11
dt2 = date(2001,1,1)
r = dt1:month(1):dt2
Base.Test.@test length(r) == 13
Base.Test.@test last(r) == date(2001,1,1)
Base.Test.@test typeof(r.step) == Month{ISOCalendar}
r = dt1:weeks(2):dt2
Base.Test.@test length(r) == 27
Base.Test.@test last(r) == date(2000,12,30)
Base.Test.@test typeof(r.step) == Week{ISOCalendar}
dt2 = date(2000,3,1)
r = dt2:day(-1):dt1
Base.Test.@test length(r) == 61
Base.Test.@test last(r) == date(2000,1,1)
Base.Test.@test typeof(r.step) == Day{ISOCalendar}
Base.Test.@test last(r + year(1)) == date(2000,12,31)
dt = date(2009,1,1)
dt2 = date(2013,1,1)
MemorialDay = x->month(x)==May && dayofweek(x)==Monday && dayofweekinmonth(x)==daysofweekinmonth(x)
f = dt:MemorialDay:dt2
Base.Test.@test typeof(f) == DateRange1{ISOCalendar}
Base.Test.@test [f] == [date(2009,5,25),date(2010,5,31),date(2011,5,30),date(2012,5,28)]
Base.Test.@test length(f) == 4
Base.Test.@test last(f) == date(2012,5,28)
Base.Test.@test typeof(f.step) == Function
Thanksgiving = x->dayofweek(x)==Thursday && month(x)==November && dayofweekinmonth(x)==4
f = dt:Thanksgiving:dt2
Base.Test.@test [f] == [date(2009,11,26),date(2010,11,25),date(2011,11,24),date(2012,11,22)]
Base.Test.@test length(f) == 4
Base.Test.@test last(f) == date(2012,11,22)
ff = recur(dt:dt2) do x
	dayofweek(x)==4 &&
	month(x)==11 &&
	dayofweekinmonth(x)==4
end
Base.Test.@test [ff] == [f]
#DateTime
dt = datetime(2013,7,1)
Base.Test.@test year(dt) == 2013
Base.Test.@test month(dt) == 7
Base.Test.@test day(dt) == 1
Base.Test.@test hour(dt) == 0
Base.Test.@test minute(dt) == 0
Base.Test.@test second(dt) == 0
Base.Test.@test string(dt) == "2013-07-01T00:00:00 UTC"
Base.Test.@test Datetime._days(dt) == 735050 == fld(int64(dt),86400000)
Base.Test.@test week(dt) == 27
Base.Test.@test isdatetime(dt)
Base.Test.@test calendar(dt) == ISOCalendar
Base.Test.@test timezone(dt) == Datetime.Zone0
Base.Test.@test typemax(DateTime) == datetime(292277024,12,31,23,59,59)
Base.Test.@test typemin(DateTime) == datetime(-292277023,1,1,0,0,0)
Base.Test.@test isleap(dt) == false
Base.Test.@test lastdayofmonth(dt) == 31
Base.Test.@test dayofweek(dt) == 1
Base.Test.@test dayofweekinmonth(dt) == 1
Base.Test.@test daysofweekinmonth(dt) == 5
Base.Test.@test firstdayofweek(dt) == dt
Base.Test.@test lastdayofweek(dt) == dt + days(6)
Base.Test.@test dayofyear(dt) == 182
Base.Test.@test typeof(now()) <: DateTime
Base.Test.@test int64(dt) == 63508320025000
Base.Test.@test second(datetime(1972,6,30,22,58,60)) == 0 #entering "invalid" periods just rolls the date forward
dt = datetime(2012,2,29)
dt2 = datetime(2000,2,1)
Base.Test.@test dt > dt2
Base.Test.@test dt != dt2
Base.Test.@test isdatetime(dt)
Base.Test.@test calendar(dt) == ISOCalendar
Base.Test.@test year(-dt) == -2012
Base.Test.@test +dt == dt
Base.Test.@test_throws dt + dt2
Base.Test.@test_throws dt * dt2
Base.Test.@test_throws dt / dt2
Base.Test.@test dt - dt2 == 381110402000
Base.Test.@test dt > dt2
Base.Test.@test dt2 < dt
Base.Test.@test dt2 - year(3) == datetime(1997,2,1)
Base.Test.@test dt + year(1) == datetime(2013,2,28)
Base.Test.@test year(3) - dt2 == dt2 - year(3)
Base.Test.@test dt2 - months(3) == datetime(1999,11,1)
Base.Test.@test months(11) + dt == datetime(2013,1,29)
Base.Test.@test months(8) + dt == dt + months(8)
Base.Test.@test dt2 + days(4411) == datetime(2012,2,29)
Base.Test.@test dt2 + days(4412) == datetime(2012,3,1)
Base.Test.@test dt2 + weeks(52) == datetime(2001,1,30)
Base.Test.@test dt2 + weeks(104) == datetime(2002,1,29)
Base.Test.@test dt2 - year(4) + days(366) == datetime(1997,2,1)
Base.Test.@test dt2 - (year(4) - days(366)) == datetime(1997,2,2) #non-associative
dt = datetime(1999,12,27,0,0,0)
check = (52,52,52,52,52,52,52,1,1,1,1,1,1,1,2,2,2,2,2,2,2)
for i = 1:21
	Base.Test.@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = datetime(2000,12,25,0,0,0)
for i = 1:21
	Base.Test.@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = datetime(2030,12,23,0,0,0)
for i = 1:21
	Base.Test.@test week(dt) == check[i]
	dt = dt + day(1)
end
dt = datetime(2004,12,20,0,0,0)
check = (52,52,52,52,52,52,52,53,53,53,53,53,53,53,1,1,1,1,1,1,1)
for i = 1:21
	Base.Test.@test week(dt) == check[i]
	dt = dt + day(1)
end

#Leap second support
# _lm =
# [1972 6 30 23 59 59;
# 1972 12 31 23 59 59;
# 1973 12 31 23 59 59;
# 1974 12 31 23 59 59;
# 1975 12 31 23 59 59;
# 1976 12 31 23 59 59;
# 1977 12 31 23 59 59;
# 1978 12 31 23 59 59;
# 1979 12 31 23 59 59;
# 1981 6 30 23 59 59;
# 1982 6 30 23 59 59;
# 1983 6 30 23 59 59;
# 1985 6 30 23 59 59;
# 1987 12 31 23 59 59;
# 1989 12 31 23 59 59;
# 1990 12 31 23 59 59;
# 1992 6 30 23 59 59;
# 1993 6 30 23 59 59;
# 1994 6 30 23 59 59;
# 1995 12 31 23 59 59;
# 1997 6 30 23 59 59;
# 1998 12 31 23 59 59;
# 2005 12 31 23 59 59;
# 2008 12 31 23 59 59;
# 2012 6 30 23 59 59]
# #Generate leap moments in microsends
# _leapmoments = ref(Int64) #1 sec before a leap second in a timeline that ignores leap seconds
# _leapmoments1 = ref(Int64) #leap second moment in leap second timeline
# for i in 1:size(_lm)[1]
# 	l = 1000*(_lm[i,6] +
# 		60*_lm[i,5] + 3600*_lm[i,4] +
# 		86400*totaldays(_lm[i,1],_lm[i,2],_lm[i,3]))
# 	push!(_leapmoments,l)
# 	push!(_leapmoments1,l+(1000*length(_leapmoments)))
# end
q = datetime(1972,6,30,23,59,58)
t = datetime(1972,6,30,23,59,59)
Base.Test.@test t - q == 1000
Base.Test.@test t == t
r = datetime(1972,7,1,0,0,0)
Base.Test.@test r - t == 2000 #Includes leap second
n = datetime(1972,6,30,23,59,60)
Base.Test.@test n - t == 1000
Base.Test.@test r - n == 1000
a = datetime(1972,12,31,23,59,59)
b = datetime(1972,12,31,23,59,60)
c = datetime(1973,1,1,0,0,0)
d = datetime(1973,1,1,0,0,1)
Base.Test.@test c - b == 1000
Base.Test.@test b - a == 1000
a = datetime(1982,6,30,23,59,59)
b = datetime(1982,6,30,23,59,60)
c = datetime(1982,7,1,0,0,0)
Base.Test.@test c - b == 1000
Base.Test.@test b - a == 1000

#Offsets
TZ = CST
a = datetime(1972,12,31,17,59,59,0,TZ)
b = datetime(1972,12,31,17,59,60,0,TZ)
c = datetime(1972,12,31,18,0,0,0,TZ)
d = datetime(1972,12,31,18,0,1,0,TZ)
Base.Test.@test c - b == 1000
Base.Test.@test b - a == 1000
Base.Test.@test c - a == 2000 #Includes leap second
bb = timezone(b,PST)
Base.Test.@test bb == b
bb = timezone(b,"Pacific/Honolulu")
Base.Test.@test bb == b
bb = offset(b,Offset{-360})
Base.Test.@test bb == b
Base.Test.@test string(bb) == "1972-12-31T17:59:60 -06:00"
aa = offset(a,Offset{-360})
Base.Test.@test aa == a
Base.Test.@test offset(a,Offset{-600}) == timezone(a,"Pacific/Honolulu")
Base.Test.@test offset(a,offset(hours(10))) == timezone(a,"Pacific/Honolulu")
Base.Test.@test timezone(datetime(2013,7,6,0,0,0,0,"America/Chicago")) == CDT

#Daylight savings time
t1 = datetime(1920,6,13,1,59,59,0,CST)
t2 = datetime(1920,6,13,2,0,0,0,CST)
t3 = datetime(1920,6,13,3,0,0,0,CST)
Base.Test.@test t2 - t1 == 1000
Base.Test.@test t3 - t1 == 1000
Base.Test.@test t3 == t2
Base.Test.@test t1 + second(1) == t2
Base.Test.@test t1 + second(1) == t3
t4 = datetime(1920,6,13,4,0,0,0,CST)
Base.Test.@test t4 - hour(1) == t3

#Conversion from Unix time
# _ = time()
# tt = unix2datetime(1000*int64(_),TZ)
# ttt = TmStruct(_)
# Base.Test.@test year(tt) == int("20"*string(digits(ttt.year)[2])*string(digits(ttt.year)[1]))
# Base.Test.@test month(tt) == ttt.month+1
# Base.Test.@test day(tt) == ttt.mday
# Base.Test.@test hour(tt) == ttt.hour
# Base.Test.@test minute(tt) == ttt.min
# Base.Test.@test second(tt) == ttt.sec #second may be off by 1, not sure why
#DateTimeRange
dt1 = datetime(2000,1,1)
dt2 = datetime(2010,1,1)
y = year(1)
r = dt1:y:dt2
Base.Test.@test first(r) == datetime(2000,1,1)
Base.Test.@test last(r) == datetime(2010,1,1)
Base.Test.@test typeof(r.step) == Year{ISOCalendar}
Base.Test.@test length(r) == 11
dt2 = datetime(2001,1,1)
r = dt1:month(1):dt2
Base.Test.@test length(r) == 13
Base.Test.@test last(r) == datetime(2001,1,1)
Base.Test.@test typeof(r.step) == Month{ISOCalendar}
r = dt1:weeks(2):dt2
Base.Test.@test length(r) == 27
Base.Test.@test last(r) == datetime(2000,12,30)
Base.Test.@test typeof(r.step) == Week{ISOCalendar}
dt2 = datetime(2000,3,1)
r = dt2:day(-1):dt1
Base.Test.@test length(r) == 61
Base.Test.@test last(r) == datetime(2000,1,1)
Base.Test.@test typeof(r.step) == Day{ISOCalendar}
Base.Test.@test last(r + year(1)) == datetime(2000,12,31)
dt2 = datetime(2000,1,31)
r = dt1:hours(24):dt2
Base.Test.@test length(r) == 31
Base.Test.@test last(r) == datetime(2000,1,31)
Base.Test.@test typeof(r.step) == Hour{ISOCalendar}
r = dt1:minutes(1440):dt2
Base.Test.@test length(r) == 31
Base.Test.@test last(r) == datetime(2000,1,31)
Base.Test.@test typeof(r.step) == Minute{ISOCalendar}
r = dt1:seconds(86400):dt2
Base.Test.@test length(r) == 31
Base.Test.@test last(r) == datetime(2000,1,31)
Base.Test.@test typeof(r.step) == Second{ISOCalendar}

ext: jl


text:
# This script parses and plots a real time series from a ROS bag where.
# Bags are loaded from a given directory and *should be paired with yaml files
# of the same name.* This script then parses the yaml files to identify ranges
# of values for given params and succinctly selects and processes sets of
# repeated trials.

using PyPlot
using RosDataProcess

# Configuration
save_output = false
colors = generate_colors(10)

println("Loading bag directory")
bag_directory = "$(homedir())/bagfiles/dynamic_exploration_journal/"
bags = load_directory(bag_directory)

# show topics
println("Topic names:")
foreach(println, get_topic_names(bags[1]))

params = Dict(:num_decentralized =>
              "dynamic_exploration/num_decentralized_planning_rounds")

ranges = Dict(key => get_range(value, bags) for (key, value) in params)

# Plot entropy results.
figure()
for (ii, num_decentralized) = enumerate(ranges[:num_decentralized])
  println("Plotting trials:"
          * " num_decentralized=$(num_decentralized)")

  println("  Filtering")
  fixed = Dict(params[:num_decentralized] => num_decentralized)
  trials = filter_bags(fixed, bags)

  println("  Reading bags")
  entropy = read_topic("/entropy", trials; accessor=x->x.data,
                       interpolate=true)

  println("  Plotting data")
  plots = plot_trials(entropy, mean=true, standard_error=true,
                      trials=false, color=colors[ii])

  legend_string(n) = (n == 0 ? "G" : "D\$_$(n)\$")
  plots[:mean][1].set_label(legend_string(num_decentralized))
end

title("Entropy reduction: dynamic exploration")
xlabel("Time (s)")
ylabel("Objective (bits)")
legend(loc="lower right")

if save_output
  save_latex("./fig", "dynamic_entropy_reduction")
end

foreach(close, bags)

max_stars_repo_name: mcorah/RosDataProcess

max_issues_repo_name: mcorah/RosDataProcess

max_forks_repo_name: mcorah/RosDataProcess


text:
module Shearheating

# This implements different methods to specify Shearheating of rocks
#
# If you want to add a new method here, feel free to do so.
# Remember to also export the function name in GeoParams.jl (in addition to here)

using Parameters, LaTeXStrings, Unitful
using ..Units
using GeoParams: AbstractMaterialParam
import Base.show

abstract type AbstractShearheating <: AbstractMaterialParam end

export  ComputeShearheating, ComputeShearheating!,  # calculation routines
        ConstantShearheating                        # constant

# Constant Shearheating -------------------------------------------------------
"""
    ConstantShearheating(Χ=0.0NoUnits)

Set the shear heating efficiency [0-1] parameter
```math
Χ  = cst
```
where ``\\Chi`` is the shear heating efficiency [NoUnits]

Shear heating is computed as
```math
H_s = \\Chi \\cdot \\tau_{ij}(\\dot{\\varepsilon}_{ij} - \\dot{\\varepsilon}^{el}_{ij})
```

"""
@with_kw_noshow mutable struct ConstantShearheating <: AbstractShearheating
    equation::LaTeXString   =   L"\H_s = \Chi \tau_{ij}(\dot{\varepsilon}_{ij} - \dot{\varepsilon}^{el}_{ij})"
    Χ::GeoUnit              =   0NoUnits
end

ConstantShearheating(x::Number) = ConstantShearheating(Χ = x*NoUnits)   # allow

# In-place routine
function ComputeShearheating!(H_s, τ, ε, ε_el, s::ConstantShearheating)
    @unpack Χ   = s

    if isnothing(ε_el)
        H_s = Value(Χ)*sum( τ .* ε  )
    else
        H_s = Value(Χ)*sum( τ .* (ε .- ε_el) )
    end

end

# Calculation routine
function ComputeShearheating(τ, ε, ε_el, s::ConstantShearheating)
    @unpack Χ   = s

    if isnothing(ε_el)
        H_s = Value(Χ)*sum( τ .* ε  )
    else
        H_s = Value(Χ)*sum( τ .* (ε .- ε_el) )
    end

    return H_s
end

# Print info
function show(io::IO, g::ConstantShearheating)
    print(io, "Shear heating: H_s = $(g.Χ.val) τ_ij*(ε_ij - ε^el_ij)")
end
#-------------------------------------------------------------------------

# Help info for the calculation routines
"""
    H_s = ComputeShearheating(τ, ε, ε_el,  s:<AbstractShearheating)

Computes the shear heating source term

```math
H_s = \\Chi \\cdot \\tau_{ij} ( \\dot{\\varepsilon}_{ij} - \\dot{\\varepsilon}^{el}_{ij})
```

# Parameters
- ``\\Chi`` : The efficiency of shear heating (between 0-1)
- ``\\tau_{ij}`` : The full deviatoric stress tensor [4 components in 2D; 9 in 3D]
- ``\\dot{\\varepsilon}_{ij}`` : The full deviatoric strainrate tensor
- ``\\dot{\\varepsilon}^{el}_{ij}`` : The full elastic deviatoric strainrate tensor

"""
ComputeShearheating(τ, ε, ε_el, s::AbstractShearheating)

"""
    H_s = ComputeShearheating(τ, ε,  s:<AbstractShearheating)

Computes the shear heating source term when there is no elasticity

```math
H_s = \\Chi \\cdot \\tau_{ij}  \\dot{\\varepsilon}_{ij}
```

# Parameters
- ``\\Chi`` : The efficiency of shear heating (between 0-1)
- ``\\tau_{ij}`` : The full deviatoric stress tensor [4 components in 2D; 9 in 3D]
- ``\\dot{\\varepsilon}_{ij}`` : The full deviatoric strainrate tensor
"""
ComputeShearheating(τ::Any, ε::Any, s::AbstractShearheating) = ComputeShearheating(τ, ε, nothing, s::ConstantShearheating)

"""
    ComputeShearheating!(H_s, τ, ε, ε_el,  s:<AbstractShearheating)

Computes the shear heating source term in-place

```math
H_s = \\Chi \\cdot \\tau_{ij} ( \\dot{\\varepsilon}_{ij} - \\dot{\\varepsilon}^{el}_{ij})
```

# Parameters
- ``\\Chi`` : The efficiency of shear heating (between 0-1)
- ``\\tau_{ij}`` : The full deviatoric stress tensor [4 components in 2D; 9 in 3D]
- ``\\dot{\\varepsilon}_{ij}`` : The full deviatoric strainrate tensor
- ``\\dot{\\varepsilon}^{el}_{ij}`` : The full elastic deviatoric strainrate tensor

*NOTE:*
The shear heating terms require the full deviatoric stress & strain rate tensors, i.e.:

```math
2D: \\tau_{ij} = \\left(
                \\begin{matrix}
                    \\tau_{xx} & \\tau_{xz} \\\\
                    \\tau_{zx} & \\tau_{zz}
                \\end{matrix}
            \\right)
```
Since ``\\tau_{zx}=\\tau_{xz}``, most geodynamic codes only take one of the terms into account; shear heating requires all components to be used!
"""
ComputeShearheating!(H_s, τ, ε, ε_el, s::AbstractShearheating)

"""
    ComputeShearheating!(H_s, τ, ε,  s:<AbstractShearheating)

Computes the shear heating source term `H_s` in-place when there is no elasticity

```math
H_s = \\Chi \\cdot \\tau_{ij}  \\dot{\\varepsilon}_{ij}
```

# Parameters
- ``\\Chi`` : The efficiency of shear heating (between 0-1)
- ``\\tau_{ij}`` : The full deviatoric stress tensor [4 components in 2D; 9 in 3D]
- ``\\dot{\\varepsilon}_{ij}`` : The full deviatoric strainrate tensor
"""
ComputeShearheating!(H_s::Any, τ::Any, ε::Any, s::AbstractShearheating) = ComputeShearheating!(H_s, τ, ε, nothing, s::ConstantShearheating)

end


text:
using Oceananigans.Architectures: device_event
using Oceananigans: short_show

struct FFTBasedPoissonSolver{A, G, Λ, S, B, T}
    architecture :: A
            grid :: G
     eigenvalues :: Λ
         storage :: S
          buffer :: B
      transforms :: T
end

transform_str(transform) = string(typeof(transform).name.wrapper, ", ")

function transform_list_str(transform_list)
    transform_strs = (transform_str(t) for t in transform_list)
    list = string(transform_strs...)
    list = list[1:end-2]
    return list
end

Base.show(io::IO, solver::FFTBasedPoissonSolver{A, G}) where {A, G}=
    print(io, "FFTBasedPoissonSolver{$A}: \n",
          "├── grid: $(short_show(solver.grid))\n",
          "├── storage: $(typeof(solver.storage))\n",
          "├── buffer: $(typeof(solver.buffer))\n",
          "└── transforms:\n",
          "    ├── forward: ", transform_list_str(solver.transforms.forward), "\n",
          "    └── backward: ", transform_list_str(solver.transforms.backward))

function FFTBasedPoissonSolver(arch, grid, planner_flag=FFTW.PATIENT)
    topo = (TX, TY, TZ) =  topology(grid)

    λx = poisson_eigenvalues(grid.Nx, grid.Lx, 1, TX())
    λy = poisson_eigenvalues(grid.Ny, grid.Ly, 2, TY())
    λz = poisson_eigenvalues(grid.Nz, grid.Lz, 3, TZ())

    eigenvalues = (
        λx = arch_array(arch, λx),
        λy = arch_array(arch, λy),
        λz = arch_array(arch, λz)
    )

    storage = arch_array(arch, zeros(complex(eltype(grid)), size(grid)...))

    transforms = plan_transforms(arch, grid, storage, planner_flag)

    # Need buffer for index permutations and transposes.
    buffer_needed = arch isa GPU && Bounded in topo
    buffer = buffer_needed ? similar(storage) : nothing

    return FFTBasedPoissonSolver(arch, grid, eigenvalues, storage, buffer, transforms)
end

"""
    solve!(ϕ, poisson_solver::FFTBasedPoissonSolver, b, m=0)

Solves the "generalized" Poisson equation,

```math
(∇² + m) ϕ = b,
```

where ``m`` is a number, using a eigenfunction expansion of the discrete Poisson operator
on a staggered grid and for periodic or Neumann boundary conditions.

In-place transforms are applied to ``b``, which means ``b`` must have complex-valued
elements (typically the same type as `poisson_solver.storage`).

Note: ``(∇² + m) ϕ = b`` is sometimes called the "screened Poisson" equation
when ``m < 0``, or the Helmholtz equation when ``m > 0``.
"""
function solve!(ϕ, solver::FFTBasedPoissonSolver, b, m=0)
    arch = solver.architecture
    topo = TX, TY, TZ = topology(solver.grid)
    Nx, Ny, Nz = size(solver.grid)
    λx, λy, λz = solver.eigenvalues

    # Temporarily store the solution in ϕc
    ϕc = solver.storage

    # Transform b *in-place* to eigenfunction space
    [transform!(b, solver.buffer) for transform! in solver.transforms.forward]

    # Solve the discrete screened Poisson equation (∇² + m) ϕ = b.
    @. ϕc = - b / (λx + λy + λz - m)

    # If m === 0, the "zeroth mode" at `i, j, k = 1, 1, 1` is undetermined;
    # we set this to zero by default. Another slant on this "problem" is that
    # λx[1, 1, 1] + λy[1, 1, 1] + λz[1, 1, 1] = 0, which yields ϕ[1, 1, 1] = Inf or NaN.
    m === 0 && CUDA.@allowscalar ϕc[1, 1, 1] = 0

    # Apply backward transforms in order
    [transform!(ϕc, solver.buffer) for transform! in solver.transforms.backward]

    copy_event = launch!(arch, solver.grid, :xyz, copy_real_component!, ϕ, ϕc, dependencies=device_event(arch))
    wait(device(arch), copy_event)

    return ϕ
end

@kernel function copy_real_component!(ϕ, ϕc)
    i, j, k = @index(Global, NTuple)
    @inbounds ϕ[i, j, k] = real(ϕc[i, j, k])
end


text:
function jt04()
#
#  test polynomial suggested by Jenkins and Traub
#
   p = poly([0.1;0.1;0.1;0.5;0.6;0.7]);
   z = [0.5  1; 0.6  1; 0.7  1; 0.1  3];
   p, PolyZeros(z)
end


text:
using Molly
using CUDA

using Base.Threads
using Statistics
using Test

@warn "This file does not include all the tests for Molly.jl due to CI time limits, " *
        "see the test directory for more"

if nthreads() > 1
    @info "The parallel tests will be run as Julia is running on $(nthreads()) threads"
else
    @warn "The parallel tests will not be run as Julia is running on 1 thread"
end

if CUDA.functional()
    @info "The GPU tests will be run as CUDA is available"
else
    @warn "The GPU tests will not be run as CUDA is not available"
end

CUDA.allowscalar(false) # Check that we never do scalar indexing on the GPU

@testset "Spatial" begin
    @test vector1D(4.0, 6.0, 10.0) ==  2.0
    @test vector1D(1.0, 9.0, 10.0) == -2.0
    @test vector1D(6.0, 4.0, 10.0) == -2.0
    @test vector1D(9.0, 1.0, 10.0) ==  2.0

    @test vector(SVector(4.0, 1.0, 6.0),
                    SVector(6.0, 9.0, 4.0), 10.0) == SVector(2.0, -2.0, -2.0)

    @test adjust_bounds(8.0 , 10.0) == 8.0
    @test adjust_bounds(12.0, 10.0) == 2.0
    @test adjust_bounds(-2.0, 10.0) == 8.0

    s = Simulation(
        simulator=VelocityVerlet(),
        atoms=[Atom(), Atom(), Atom()],
        coords=[SVector(1.0, 1.0, 1.0), SVector(2.0, 2.0, 2.0),
                SVector(5.0, 5.0, 5.0)],
        box_size=10.0,
        neighbour_finder=DistanceNeighbourFinder(trues(3, 3), 10, 2.0)
    )
    find_neighbours!(s, s.neighbour_finder, 0, parallel=false)
    @test s.neighbours == [(2, 1)]
    if nthreads() > 1
        find_neighbours!(s, s.neighbour_finder, 0, parallel=true)
        @test s.neighbours == [(2, 1)]
    end
end

temp = 298
timestep = 0.002
n_steps = 20_000
box_size = 2.0

@testset "Lennard-Jones gas 2D" begin
    n_atoms = 20

    s = Simulation(
        simulator=VelocityVerlet(),
        atoms=[Atom(attype="Ar", name="Ar", resnum=i, resname="Ar", charge=0.0,
                    mass=10.0, σ=0.3, ϵ=0.2) for i in 1:n_atoms],
        general_inters=(LennardJones(true),),
        coords=[box_size .* rand(SVector{2}) for i in 1:n_atoms],
        velocities=[velocity(10.0, temp, dims=2) .* 0.01 for i in 1:n_atoms],
        temperature=temp,
        box_size=box_size,
        neighbour_finder=DistanceNeighbourFinder(trues(n_atoms, n_atoms), 10, 2.0),
        thermostat=AndersenThermostat(10.0),
        loggers=Dict("temp" => TemperatureLogger(100),
                     "coords" => CoordinateLogger(100, dims=2)),
        timestep=timestep,
        n_steps=n_steps
    )

    show(devnull, s)

    @time simulate!(s, parallel=false)
end

@testset "Lennard-Jones gas" begin
    n_atoms = 100
    parallel_list = nthreads() > 1 ? (false, true) : (false,)

    for parallel in parallel_list
        s = Simulation(
            simulator=VelocityVerlet(),
            atoms=[Atom(attype="Ar", name="Ar", resnum=i, resname="Ar", charge=0.0,
                        mass=10.0, σ=0.3, ϵ=0.2) for i in 1:n_atoms],
            general_inters=(LennardJones(true),),
            coords=[box_size .* rand(SVector{3}) for i in 1:n_atoms],
            velocities=[velocity(10.0, temp) .* 0.01 for i in 1:n_atoms],
            temperature=temp,
            box_size=box_size,
            neighbour_finder=DistanceNeighbourFinder(trues(n_atoms, n_atoms), 10, 2.0),
            thermostat=AndersenThermostat(10.0),
            loggers=Dict("temp" => TemperatureLogger(100),
                         "coords" => CoordinateLogger(100),
                         "energy" => EnergyLogger(100)),
            timestep=timestep,
            n_steps=n_steps
        )

        @time simulate!(s, parallel=parallel)

        final_coords = last(s.loggers["coords"].coords)
        @test minimum(minimum.(final_coords)) > 0.0
        @test maximum(maximum.(final_coords)) < box_size
        displacements(final_coords, box_size)
        distances(final_coords, box_size)
        rdf(final_coords, box_size)
    end
end

@testset "Lennard-Jones gas velocity-free" begin
    n_atoms = 100
    coords = [box_size .* rand(SVector{3}) for i in 1:n_atoms]

    s = Simulation(
        simulator=VelocityFreeVerlet(),
        atoms=[Atom(attype="Ar", name="Ar", resnum=i, resname="Ar", charge=0.0,
                    mass=10.0, σ=0.3, ϵ=0.2) for i in 1:n_atoms],
        general_inters=(LennardJones(true),),
        coords=coords,
        velocities=[c .+ 0.01 .* rand(SVector{3}) for c in coords],
        temperature=temp,
        box_size=box_size,
        neighbour_finder=DistanceNeighbourFinder(trues(n_atoms, n_atoms), 10, 2.0),
        thermostat=AndersenThermostat(10.0),
        loggers=Dict("temp" => TemperatureLogger(100),
                     "coords" => CoordinateLogger(100)),
        timestep=timestep,
        n_steps=n_steps
    )

    @time simulate!(s, parallel=false)
end

@testset "Diatomic molecules" begin
    n_atoms = 100
    coords = [box_size .* rand(SVector{3}) for i in 1:(n_atoms / 2)]
    for i in 1:length(coords)
        push!(coords, coords[i] .+ [0.1, 0.0, 0.0])
    end
    bonds = [HarmonicBond(i, i + (n_atoms ÷ 2), 0.1, 300_000.0) for i in 1:(n_atoms ÷ 2)]

    s = Simulation(
        simulator=VelocityVerlet(),
        atoms=[Atom(attype="H", name="H", resnum=i, resname="H", charge=0.0,
                    mass=10.0, σ=0.3, ϵ=0.2) for i in 1:n_atoms],
        specific_inter_lists=(bonds,),
        general_inters=(LennardJones(true),),
        coords=coords,
        velocities=[velocity(10.0, temp) .* 0.01 for i in 1:n_atoms],
        temperature=temp,
        box_size=box_size,
        neighbour_finder=DistanceNeighbourFinder(trues(n_atoms, n_atoms), 10, 2.0),
        thermostat=AndersenThermostat(10.0),
        loggers=Dict("temp" => TemperatureLogger(10),
                        "coords" => CoordinateLogger(10)),
        timestep=timestep,
        n_steps=n_steps
    )

    @time simulate!(s, parallel=false)
end

@testset "Peptide" begin
    timestep = 0.0002
    n_steps = 100
    atoms, specific_inter_lists, general_inters, nb_matrix, coords, box_size = readinputs(
                normpath(@__DIR__, "..", "data", "5XER", "gmx_top_ff.top"),
                normpath(@__DIR__, "..", "data", "5XER", "gmx_coords.gro"))

    true_n_atoms = 5191
    @test length(atoms) == true_n_atoms
    @test length(coords) == true_n_atoms
    @test size(nb_matrix) == (true_n_atoms, true_n_atoms)
    @test length(specific_inter_lists) == 3
    @test length(general_inters) == 2
    @test box_size == 3.7146
    show(devnull, first(atoms))

    s = Simulation(
        simulator=VelocityVerlet(),
        atoms=atoms,
        specific_inter_lists=specific_inter_lists,
        general_inters=general_inters,
        coords=coords,
        velocities=[velocity(a.mass, temp) .* 0.01 for a in atoms],
        temperature=temp,
        box_size=box_size,
        neighbour_finder=DistanceNeighbourFinder(nb_matrix, 10),
        thermostat=AndersenThermostat(10.0),
        loggers=Dict("temp" => TemperatureLogger(10),
                        "coords" => CoordinateLogger(10),
                        "energy" => EnergyLogger(10)),
        timestep=timestep,
        n_steps=n_steps
    )

    @time simulate!(s, parallel=false)
end

@testset "Float32" begin
    timestep = 0.0002f0
    n_steps = 100
    atoms, specific_inter_lists, general_inters, nb_matrix, coords, box_size = readinputs(
                Float32,
                normpath(@__DIR__, "..", "data", "5XER", "gmx_top_ff.top"),
                normpath(@__DIR__, "..", "data", "5XER", "gmx_coords.gro"))

    s = Simulation(
        simulator=VelocityVerlet(),
        atoms=atoms,
        specific_inter_lists=specific_inter_lists,
        general_inters=general_inters,
        coords=coords,
        velocities=[velocity(Float32, a.mass, temp) .* 0.01f0 for a in atoms],
        temperature=Float32(temp),
        box_size=box_size,
        neighbour_finder=DistanceNeighbourFinder(nb_matrix, 10, 1.2f0),
        thermostat=AndersenThermostat(10.0f0),
        loggers=Dict("temp" => TemperatureLogger(Float32, 10),
                        "coords" => CoordinateLogger(Float32, 10),
                        "energy" => EnergyLogger(Float32, 10)),
        timestep=timestep,
        n_steps=n_steps
    )

    @time simulate!(s, parallel=false)
end

@testset "General interactions" begin
    n_atoms = 100
    general_inter_types = (
        LennardJones(true ), SoftSphere(true ), Mie(5, 10, true ), Coulomb(true ), Gravity(10.0, true ),
        LennardJones(false), SoftSphere(false), Mie(5, 10, false), Coulomb(false), Gravity(10.0, false),
    )

    @testset "$gi" for gi in general_inter_types
        if gi.nl_only
            neighbour_finder = DistanceNeighbourFinder(trues(n_atoms, n_atoms), 10, 1.5)
        else
            neighbour_finder = NoNeighbourFinder()
        end

        s = Simulation(
            simulator=VelocityVerlet(),
            atoms=[Atom(charge=i % 2 == 0 ? -1.0 : 1.0, mass=10.0, σ=0.2, ϵ=0.2) for i in 1:n_atoms],
            general_inters=(gi,),
            coords=placeatoms(n_atoms, box_size, 0.2),
            velocities=[velocity(10.0, temp) .* 0.01 for i in 1:n_atoms],
            temperature=temp,
            box_size=box_size,
            neighbour_finder=neighbour_finder,
            thermostat=AndersenThermostat(10.0),
            loggers=Dict("temp" => TemperatureLogger(100),
                         "coords" => CoordinateLogger(100),
                         "energy" => EnergyLogger(100)),
            timestep=timestep,
            n_steps=n_steps
        )

        @time simulate!(s)
    end
end

@enum Status susceptible infected recovered

# Custom atom type
mutable struct Person
    status::Status
    mass::Float64
    σ::Float64
    ϵ::Float64
end

# Custom GeneralInteraction
struct SIRInteraction <: GeneralInteraction
    nl_only::Bool
    dist_infection::Float64
    prob_infection::Float64
    prob_recovery::Float64
end

# Custom Logger
struct SIRLogger <: Logger
    n_steps::Int
    fracs_sir::Vector{Vector{Float64}}
end

@testset "Different implementations" begin
    function placediatomics(n_molecules, box_size, min_dist, bond_length)
        min_dist_sq = min_dist ^ 2
        coords = SArray{Tuple{3}, Float64, 1, 3}[]
        while length(coords) < (n_molecules * 2)
            new_coord_a = rand(SVector{3}) .* box_size
            new_coord_b = copy(new_coord_a) + SVector{3}([bond_length, 0.0, 0.0])
            okay = new_coord_b[1] <= box_size
            for coord in coords
                if sum(abs2, vector(coord, new_coord_a, box_size)) < min_dist_sq ||
                        sum(abs2, vector(coord, new_coord_b, box_size)) < min_dist_sq
                    okay = false
                    break
                end
            end
            if okay
                push!(coords, new_coord_a)
                push!(coords, new_coord_b)
            end
        end
        return coords
    end

    n_atoms = 400
    mass = 10.0
    box_size = 6.0
    temp = 1.0
    starting_coords = placediatomics(n_atoms ÷ 2, box_size, 0.2, 0.2)
    starting_velocities = [velocity(mass, temp) for i in 1:n_atoms]
    starting_coords_f32 = [Float32.(c) for c in starting_coords]
    starting_velocities_f32 = [Float32.(c) for c in starting_velocities]

    function runsim(nl::Bool, parallel::Bool, gpu_diff_safe::Bool, f32::Bool, gpu::Bool)
        n_atoms = 400
        n_steps = 200
        mass = f32 ? 10.0f0 : 10.0
        box_size = f32 ? 6.0f0 : 6.0
        timestep = f32 ? 0.02f0 : 0.02
        temp = f32 ? 1.0f0 : 1.0
        simulator = VelocityVerlet()
        thermostat = NoThermostat()
        bonds = [HarmonicBond((i * 2) - 1, i * 2, f32 ? 0.2f0 : 0.2, f32 ? 10_000.0f0 : 10_000.0) for i in 1:(n_atoms ÷ 2)]
        specific_inter_lists = (bonds,)

        neighbour_finder = NoNeighbourFinder()
        general_inters = (LennardJones(false),)
        if nl
            neighbour_finder = DistanceNeighbourFinder(trues(n_atoms, n_atoms), 10, f32 ? 1.5f0 : 1.5)
            general_inters = (LennardJones(true),)
        end

        if gpu
            coords = cu(deepcopy(f32 ? starting_coords_f32 : starting_coords))
            velocities = cu(deepcopy(f32 ? starting_velocities_f32 : starting_velocities))
            atoms = cu([AtomMin(f32 ? 0.0f0 : 0.0, mass, f32 ? 0.2f0 : 0.2, f32 ? 0.2f0 : 0.2) for i in 1:n_atoms])
        else
            coords = deepcopy(f32 ? starting_coords_f32 : starting_coords)
            velocities = deepcopy(f32 ? starting_velocities_f32 : starting_velocities)
            atoms = [Atom(attype="Ar", name="Ar", resnum=i, resname="Ar", charge=f32 ? 0.0f0 : 0.0,
                            mass=mass, σ=f32 ? 0.2f0 : 0.2, ϵ=f32 ? 0.2f0 : 0.2) for i in 1:n_atoms]
        end

        s = Simulation(
            simulator=simulator,
            atoms=atoms,
            specific_inter_lists=specific_inter_lists,
            general_inters=general_inters,
            coords=coords,
            velocities=velocities,
            temperature=temp,
            box_size=box_size,
            neighbour_finder=neighbour_finder,
            thermostat=thermostat,
            timestep=timestep,
            n_steps=n_steps,
            gpu_diff_safe=gpu_diff_safe
        )

        c = simulate!(s, parallel=parallel)
        return c
    end

    runs = [
        ("in-place"        , [false, false, false, false, false]),
        ("in-place NL"     , [true , false, false, false, false]),
        ("in-place f32"    , [false, false, false, true , false]),
        ("out-of-place"    , [false, false, true , false, false]),
        ("out-of-place f32", [false, false, true , true , false]),
    ]
    if nthreads() > 1
        push!(runs, ("in-place parallel"   , [false, true , false, false, false]))
        push!(runs, ("in-place NL parallel", [true , true , false, false, false]))
    end
    if CUDA.functional()
        push!(runs, ("out-of-place gpu"    , [false, false, true , false, true ]))
        push!(runs, ("out-of-place gpu f32", [false, false, true , true , true ]))
    end

    final_coords_ref = Array(runsim(runs[1][2]...))
    for (name, args) in runs
        final_coords = Array(runsim(args...))
        final_coords_f64 = [Float64.(c) for c in final_coords]
        diff = sum(sum(map(x -> abs.(x), final_coords_f64 .- final_coords_ref))) / (3 * n_atoms)
        # Check all simulations give the same result to within some error
        @info "$(rpad(name, 20)) - difference per coordinate $diff"
        @test diff < 1e-4
    end
end

@testset "Agent-based modelling" begin
    # Custom force function
    function Molly.force!(forces, inter::SIRInteraction, s::Simulation, i::Integer, j::Integer)
        if i == j
            # Recover randomly
            if s.atoms[i].status == infected && rand() < inter.prob_recovery
                s.atoms[i].status = recovered
            end
        elseif (s.atoms[i].status == infected && s.atoms[j].status == susceptible) ||
                    (s.atoms[i].status == susceptible && s.atoms[j].status == infected)
            # Infect close people randomly
            dr = vector(s.coords[i], s.coords[j], s.box_size)
            r2 = sum(abs2, dr)
            if r2 < inter.dist_infection ^ 2 && rand() < inter.prob_infection
                s.atoms[i].status = infected
                s.atoms[j].status = infected
            end
        end
        return nothing
    end

    # Custom logging function
    function Molly.log_property!(logger::SIRLogger, s::Simulation, step_n::Integer)
        if step_n % logger.n_steps == 0
            counts_sir = [
                count(p -> p.status == susceptible, s.atoms),
                count(p -> p.status == infected   , s.atoms),
                count(p -> p.status == recovered  , s.atoms)
            ]
            push!(logger.fracs_sir, counts_sir ./ length(s.atoms))
        end
    end

    temp = 0.01
    timestep = 0.02
    box_size = 10.0
    n_steps = 1_000
    n_people = 500
    n_starting = 2
    atoms = [Person(i <= n_starting ? infected : susceptible, 1.0, 0.1, 0.02) for i in 1:n_people]
    coords = [box_size .* rand(SVector{2}) for i in 1:n_people]
    velocities = [velocity(1.0, temp, dims=2) for i in 1:n_people]
    general_inters = (LennardJones = LennardJones(true),
                            SIR = SIRInteraction(false, 0.5, 0.06, 0.01))

    s = Simulation(
        simulator=VelocityVerlet(),
        atoms=atoms,
        general_inters=general_inters,
        coords=coords,
        velocities=velocities,
        temperature=temp,
        box_size=box_size,
        neighbour_finder=DistanceNeighbourFinder(trues(n_people, n_people), 10, 2.0),
        thermostat=AndersenThermostat(5.0),
        loggers=Dict("coords" => CoordinateLogger(10, dims=2),
                        "SIR" => SIRLogger(10, [])),
        timestep=timestep,
        n_steps=n_steps
    )

    @time simulate!(s, parallel=false)
end


text:
module ParaDwell

using Shapefile, GeoInterface, PlotlyJS, CSV, DataFrames
using LightXML, Interact, Blink, Images, CSVFiles, StatsBase

include("structs.jl")
include("Init.jl")
include("Mapping/Coastline.jl")
include("Mapping/Conversion.jl")
include("Mapping/SummariseOsData.jl")
include("XmlTools/OrdnanceSurvey.jl")
include("XmlTools/structs.jl")
include("XmlTools/xmlParse.jl")
# include("UI/MainUI_Launch.jl")
# include("UI/ClickTile.jl")

end


text:
module Qakbot
using CRC32c,Dates

mutable struct OBJ
    index
    mt
end

function date2seed(date, seed)
    dx = trunc(Int64, (Dates.day(date)-1) / 10)
    dx = dx<=2 ? dx : 2
    data = "$dx.$(lowercase(Dates.monthabbr(date))).$(Dates.year(date)).$seed"
    crc = crc32c(data)
    return crc
end

function _int32(x)
    return UInt32(0xFFFFFFFF&x)
end

function rand_int!(obj, lower, upper)
    r = extract_num!(obj)
    r &= 0xFFFFFFF
    t = lower + r/(2^28)*(upper-lower+1)
    return trunc(Int64, t)
end

function extract_num!(obj)
    if obj.index >=625
        twist!(obj)
    end
    y = obj.mt[obj.index]
    y = y ⊻ y >> 11
    y = y ⊻ y << 7 & 2636928640
    y = y ⊻ y << 15 & 4022730752
    y = y ⊻ y >> 18

    obj.index += + 1
    return _int32(y)
end

function twist!(obj)
    for i in 1:624
        y = _int32((obj.mt[i] & 0x80000000) + (obj.mt[(i + 1)%624 + 1] & 0x7fffffff))
        obj.mt[i] = obj.mt[(i + 397)%624 + 1] ⊻ y >> 1

        if y % 2 != 0
            obj.mt[i] = obj.mt[i] ⊻ 0x9908b0df
        end
    end
    obj.index = 1
end

function gen_one(date=missing, seed=rand(Int16), tlds=missing, sandbox=false)
    domain = ""
    if date === missing
        date = Dates.Date(Dates.now())
    end
    obj = OBJ(624, zero(1:624))
    obj.mt[1] = date2seed(date, seed) + sandbox
    for i in 2:624
        obj.mt[i] = _int32(1812433253 * (obj.mt[i-1] ⊻ obj.mt[i-1]>>30) + i-1)
    end
    len = rand_int!(obj,8,25)
    for l in 1:len
        domain *= 'a'+rand_int!(obj,0,25)
    end
    if tlds !== missing && length(tlds) != 0
        tld_nr = rand_int!(obj,0,length(tlds)-1)
        domain *= "." * tlds[tld_nr+1]
    end

    return domain
end

function gen_many(cnt=100)
    res = []
    date = Dates.Date(Dates.now())
    seed = rand(typeof(cnt))
    obj = OBJ(624, zero(1:624))
    obj.mt[1] = date2seed(date, seed)
    for i in 2:624
        obj.mt[i] = _int32(1812433253 * (obj.mt[i-1] ⊻ obj.mt[i-1]>>30) + i-1)
    end
    for _ in 1:cnt
        domain = ""
        len = rand_int!(obj,8,25)
        for l in 1:len
            domain *= 'a'+rand_int!(obj,0,25)
        end
        res = [res..., domain]
    end

    return res
end

end  # module Qakbot


text:
const millisInOneSecond   = Int(1e3)
const microsInOnSecond    = Int(1e6)
const nanosInOneSecond    = Int(1e9)
const nanosInOneMilli     = Int(1e6)
const nanosInOneMicro     = Int(1e3)
const secondsInOneMinute  = Int(60)
const secondsInOneHour    = Int(3600)
const secondsInOneDay     = secondsInOneHour * Int(24)
const nanosInOneMinute    = nanosInOneSecond * secondsInOneMinute
const nanosInOneHour      = nanosInOneSecond * secondsInOneHour
const nanosInOneDay       = nanosInOneSecond * secondsInOneDay

const _epoch_millis = UInt64(DateTime(1970,1,1).instant.periods.value)
const _est5edt = TimeZone("EST5EDT")

function currentTimeOfDay()
  n = now()
  h = Dates.hour(n)
  m = Dates.minute(n)
  s = Dates.second(n)
  TimeOfDay(h * nanosInOneHour + m * nanosInOneMinute + s * nanosInOneSecond)
end

calcSecondsSinceMidnight(hour::Int, minute::Int) = hour * secondsInOneHour + secondsInOneMinute

# NOTE: DateTime module does NOT have timezone info, so all those calls are UTC

function getHoursAdjustmentFromUTC(year::Integer, month::Integer, day::Integer; timezone = _est5edt)
  # create a historically-aware timezone-aware datetime
  dt = ZonedDateTime(DateTime(year,month,day), timezone)

  # extract the hours offset from UTC
  round(Int, -TimeZones.value(dt.zone.offset) / 3600)
end

createUTCDateTimeFromSecondsSinceEpoch(secondsSinceEpoch::Integer) = DateTime(Dates.UTM(secondsSinceEpoch * millisInOneSecond + _epoch_millis))

function calcSecondsEpochToMidnight(secondsSinceEpoch::Integer)

  dt = createUTCDateTimeFromSecondsSinceEpoch(secondsSinceEpoch)

  # get the hour adjustment
  y = Dates.year(dt)
  m = Dates.month(dt)
  d = Dates.day(dt)
  hourAdjustment = getHoursAdjustmentFromUTC(y, m, d)

  millisMidnightUTC = UInt64(DateTime(y, m, d).instant.periods.value)
  millisMidnightEST = UInt64(millisMidnightUTC + hourAdjustment * secondsInOneHour * millisInOneSecond)

  return round(UInt64, (millisMidnightEST - _epoch_millis) / millisInOneSecond)
end

function convertSecondsSinceEpochToSecondsSinceMidnight(secondsSinceEpoch::Integer)
  secondsSinceEpoch - calcSecondsEpochToMidnight(secondsSinceEpoch)
end

#############################################
# TimeOfDay
#############################################

type NegativeTimeOfDayError{T} <: Exception
  val::T
end

immutable TimeOfDay <: Number
  nanosSinceMidnight::Int

  function TimeOfDay{T<:Number}(val::T)
    if (val < zero(T))
      throw(NegativeTimeOfDayError{T}(val))
      # error("TimeOfDay constructor is negative! $val")
    end
    new(round(Int, val))
  end
end

function TimeOfDay(str::AbstractString)
  tmp = split(str,":")
  minutes = 0
  seconds = 0
  nanos = 0

  hours = parse(Int, tmp[1])
  if length(tmp) > 1
    minutes = parse(Int, tmp[2])
  end

  if length(tmp) > 2
    tmp2 = split(tmp[3], ".")
    seconds = parse(Int, tmp2[1])
    if length(tmp2) > 1
      nanos = parse(Int, rpad(tmp2[2],9,"0")[1:9])
    end
  end

  TimeOfDay(hours * nanosInOneHour + minutes * nanosInOneMinute + seconds * nanosInOneSecond + nanos)
end

Base.convert(::Type{TimeOfDay}, str::AbstractString) = TimeOfDay(str)

Base.typemin(::Type{TimeOfDay}) = TimeOfDay(0)
Base.typemax(::Type{TimeOfDay}) = TimeOfDay(nanosInOneDay)

function hmsu(t::TimeOfDay)
    secsSinceMidnight, nanos = divrem(t.nanosSinceMidnight, nanosInOneSecond)
    hours, hourrem = divrem(secsSinceMidnight, secondsInOneHour)
    minutes, seconds = divrem(hourrem, secondsInOneMinute)
    microseconds = div(nanos, millisInOneSecond)
    hours, minutes, seconds, microseconds
end

function Base.show(io::IO, t::TimeOfDay)
    # hours, minutes, seconds, microseconds = hmsu(t)
    @printf(io,"%d:%02d:%02d.%06d",hmsu(t)...)
    # print(io, lpad(hours,2,"0"), ":", lpad(minutes,2,"0"), ":", lpad(seconds,2,"0"), ".", lpad(microseconds,6,"0"))
end
# @createIOMethods TimeOfDay

# define how numbers get promoted to TimeOfDay
Base.promote_rule{T<:Integer}(::Type{TimeOfDay}, ::Type{T}) = TimeOfDay
Base.promote_rule{T<:AbstractFloat}(::Type{TimeOfDay}, ::Type{T}) = TimeOfDay
Base.convert{T<:Integer}(::Type{TimeOfDay}, x::T) = TimeOfDay(x)
Base.convert{T<:AbstractFloat}(::Type{TimeOfDay}, x::T) = TimeOfDay(x)
Base.convert{T<:Number}(::Type{T}, x::TimeOfDay) = T(x.nanosSinceMidnight)

# now for each op, we define the op between TimeOfDay's, then the generic "promote" version of the op will handle conversion

# for op in (:(Base.typemax), :(Base.typemin))
#   @eval $op(::Type{TimeOfDay}) = TimeOfDay($op(Int64))
# end
for op in (:<, :>, :(==), :<=, :>=, :(Base.isless))
  @eval $op(t1::TimeOfDay, t2::TimeOfDay) = $op(t1.nanosSinceMidnight, t2.nanosSinceMidnight)
end
for op in (:+, :-, :(Base.rem))
  @eval $op(t1::TimeOfDay, t2::TimeOfDay) = TimeOfDay($op(t1.nanosSinceMidnight, t2.nanosSinceMidnight))
end

# some helpers for creating ranges of the form "timestart:timestep:timestop"
Base.colon(start::TimeOfDay, step::TimeOfDay, stop::TimeOfDay) = map(TimeOfDay, start.nanosSinceMidnight : step.nanosSinceMidnight : stop.nanosSinceMidnight)
Base.colon{T<:Number}(start::TimeOfDay, step::T, stop::TimeOfDay) = colon(start, TimeOfDay(step), stop)

iszero(t::TimeOfDay) = t.nanosSinceMidnight == 0

function timeofday_formatter(hrs::Number)
    h,m,s,u = hmsu(TimeOfDay(hrs * nanosInOneHour))
    @sprintf("%d:%02d:%02d",h,m,s)
end

hours(t::TimeOfDay) = Float64(t) / nanosInOneHour

@recipe f{T<:AbstractVector{TimeOfDay}}(::Type{T}, a::T) = (xformatter := timeofday_formatter; map(hours, a))
@recipe f{T<:AbstractMatrix{TimeOfDay}}(::Type{T}, a::T) = (xformatter := timeofday_formatter; map(hours, a))


text:
"""
    nodes(network::MangalNetwork, query::Pair...)

Returns the nodes that are part of a `MangalNetwork`, with an additional optional query.
"""
function nodes(network::MangalNetwork, query::Pair...)
    return nodes(Pair("network_id", network.id), query...)
end

"""
    nodes(taxon::MangalReferenceTaxon, query::Pair...)

Returns the nodes that are instance of a `MangalReferenceTaxon`, with an additional query.
"""
function nodes(taxon::MangalReferenceTaxon, query::Pair...)
    return nodes(Pair("taxonomy_id", taxon.id), query...)
end


text:
function parse_graph_data(name::String)
    local n
    E = Tuple{Int64, Int64}[]
    C = Float64[]
    open(name) do file
        while ! eof(file)
            line = readline(file)
            designator = line[begin]
            if (designator == 'p')
                elems = split(line, " ", keepempty = false)

                @assert(length(elems) == 4)

                n = parse(Int, elems[3])
            elseif (designator == 'a')
                elems = split(line, " ", keepempty = false)

                @assert(length(elems) == 4)

                v₁ = parse(Int64, elems[2])
                v₂ = parse(Int64, elems[3])
                push!(E, minmax(v₁, v₂))

                Cᵢ = parse(Int64, elems[4])
                push!(C, Cᵢ)
            end
        end
    end
    return n, E, C
end

function generate_uncertain_costs(dist::Distribution, c::Vector{Float64})
    return map(cᵢ -> cᵢ * rand(dist), c)
end

function generate_scenario(dist::Distribution, A::Vector{InputEdge})
    return [e.c + rand(dist) * e.d for e ∈ A]
end

function calculate_cost(A::Vector{InputEdge}, x::Vector{Tuple{Int64, Int64}}, S::Vector{Float64})
    @assert length(A) == length(S)
    return sum([e.C + S[i] for (i, e) ∈ enumerate(A) if (e.i, e.j) in x || (e.j, e.i) in x])
end

function calculate_cost(A::Vector{InputEdge}, x::Vector{Tuple{Int64, Int64}}, c::Float64)
    C = sum([e.C for e ∈ A if (e.i, e.j) in x || (e.j, e.i) in x]) # first stage costs
    return C + c # complete cost
end

function get_first_stage_solution(x::Array{Float64}, E::Vector{InputEdge})::Vector{Tuple{Int, Int}}
    return [(e.i, e.j) for (i, e) ∈ enumerate(E) if x[i] > 0]
end


text:
include("coloring.jl")
include("knapsack.jl")
include("tsptw.jl")
include("nqueens.jl")
# TODO: create a latin_square testset
# inlcude("latin.jl")


text:
# License for this file: MIT (expat)
# Copyright 2021, DLR Institute of System Dynamics and Control

module SilentNoPlot

include("plot.jl")

plot(result, names::AbstractMatrix; heading::AbstractString="", grid::Bool=true, xAxis="time",
     figure::Int=1, prefix::AbstractString="", reuse::Bool=false, maxLegend::Integer=10,
     minXaxisTickLabels::Bool=false, MonteCarloAsArea=false) = nothing
showFigure(figure::Int)  = nothing
closeFigure(figure::Int) = nothing
saveFigure(figure::Int, fileName::AbstractString) = nothing

"""
    closeAllFigures()

Close all figures.
"""
closeAllFigures() = nothing

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
function LC(p, var)
    d = maxdegree(p, var)
    d == 0 && return zero(p)
    return polynomial([subs(t, var => 1) for t in terms(p) if degree(t, var) == d])
end

function pseudorem(f, g, y)
    r = f
    lc = LC(g, y)
    m = maxdegree(g, y)
    d = maxdegree(r, y)

    while d ≥ m && !iszero(r)
        r = lc * r - LC(r, y) * g * y ^ (d - m)
        d = maxdegree(r, y)
    end

    return r
end

function triangulize(H, vars)
    n = length(H)

    T = similar(H)
    @inbounds for i = n:-1:1
        degrees = Int[] # degrees of polynomials in C
        indeces = Int[] # indeces in H of polynomials in C

        # initialize
        @inbounds for (idx, p) in enumerate(H)
            deg = maxdegree(p, vars[i])
            if deg > 0
                push!(degrees, deg)
                push!(indeces, idx)
            end
        end

        # If there are several polynomials in C but none has degree 1 in x, transform the
        # set until there is a polynomial of degree 1 or only one polynomial is left
        C = H[indeces]
        while minimum(degrees) > 1 && length(C) > 1
            _, idx_max = findmax(degrees)
            _, idx_min = findmin(degrees)

            # special case when all  have same degree
            if idx_max == idx_min
                idx_min = 1
                idx_max = 2
            end

            r = pseudorem(C[idx_max], C[idx_min], vars[i])
            d = maxdegree(r, vars[i])
            if iszero(d)
                H[indeces[idx_max]] = r
                deleteat!(C, idx_max)
                deleteat!(degrees, idx_max)
                deleteat!(indeces, idx_max)
            else
                C[idx_max] = r
                degrees[idx_max] = d
            end
        end

        # if only one polynomial in C, add to triangular set and go to the next variable
        if length(indeces) == 1
            T[i] = H[indeces[1]]
            deleteat!(H, indeces[1])
            continue
        end

        # if multiple polynomials but at least one has degree 1, add this to triangular set
        # and replace the others in H with the pseudoremainder
        idx = findlast(==(1), degrees)
        if !isnothing(idx)
            pivot = H[indeces[idx]]
            T[i] = pivot
            #pushfirst!(T, pivot)

            @inbounds for j in indeces
                if j != indeces[idx]
                    H[j] = pseudorem(H[j], pivot, vars[i])
                end
            end
            deleteat!(H, indeces[idx])
            continue
        end

    end
    return T
end


text:
@testset "WorldClim BioClim" begin
    using SimpleSDMDataSources: rastername, zipurl, zipname, zippath

    zip_file = "https://biogeo.ucdavis.edu/data/worldclim/v2.1/base/wc2.1_10m_bio.zip"
    @test zipurl(WorldClim{BioClim}, 2, 10.0) == zip_file
    @test zipname(WorldClim{BioClim}, 2, 10.0) == "wc2.1_10m_bio.zip"

    raster_file = joinpath(ENV["ECODATASOURCES_PATH"], "WorldClim/BioClim/wc2.1_10m_bio_2.tif")
    @test rasterpath(WorldClim{BioClim}, 2, 10.0) == raster_file
    @test download_raster(WorldClim{BioClim}; layer=2, resolution=10.0) == raster_file
    @test isfile(raster_file)
end


text:
using Glimpse
const Gl = Glimpse
Entity(dio, Gl.Spatial(), Gl.Shape(), Gl.Material(),
       Gl.ProgramTag{Gl.InstancedPeelingProgram}(),
       Gl.UniformColor(RGBAf0(1.0, 0.0, 0.0, 0.7)), Gl.Selectable(),
       Gl.PolygonGeometry(Gl.Sphere(zero(Point3f0), 1.0f0)))


text:
#!/usr/bin/env julia

# Author: Anand S Bisen <Anand.Bisen@SAS.COM>
# Description
#
# The script compiles Python 3.7+ from source to a self contained directory.
# Optionally it also installs python packages specified in `requirements.txt`
#
# Usage: build-python.jl -v <python3-version> -d <destination-directory> -r <path-to-requirements.txt>

using ArgParse
import Downloads: download

struct Arguments
        uri::AbstractString
        version::AbstractString
        requirements::Union{Nothing, AbstractString}
        destination::AbstractString
end

function Arguments(version::AbstractString,
                   requirements::Union{Nothing, AbstractString},
                   destination::AbstractString
                   )

        uri="https://www.python.org/ftp/python/$(version)/Python-$(version).tgz"
        Arguments(uri, version, requirements, destination)
end

filename(o::Arguments)=basename(o.uri)
basefilename(o::Arguments)=splitext(filename(o))[1]

function parsecmd()
        s = ArgParseSettings()
        @add_arg_table! s begin

                "--version", "-v"
                help = "python version"
                arg_type = String
                default = "3.9.6"

                "--requirements", "-r"
                help = "path to requirements.txt file"
                arg_type = String

                "--build-dir", "-d"
                help = "build directory"
                arg_type = String
                default = "/python"

        end

        r=parse_args(s)
        return Arguments(r["version"], abspath(r["requirements"]), r["build-dir"])
end

function download_python(o::Arguments)
        outpath = joinpath("/tmp", filename(o))

	@info("Downloading Python $(o.uri)")
	try
		res = download(o.uri, outpath, verbose=false)
	catch e
		@error("Error downloading $(o.uri)")
	end
end

function extract_python(o::Arguments)
	@info("Extracting Python to /tmp/$(basefilename(o))")
        cmd=`tar -zxvf /tmp/$(filename(o)) -C /tmp`
        ret=run(cmd)
end

function build_python(o::Arguments)
        cd("/tmp/$(basefilename(o))")

        # Configure
        cmd=`./configure --enable-optimizations --prefix=$(o.destination)`
        @info("Executing $(cmd)")
        r=run(cmd)
        if r.exitcode != 0
                println("Exit Code: $(r.exitcode)")
                println("Error: $(r.err)")
                @error("configure failed $(cmd)")
        end

        # build
        cmd=`make -j2`
        @info("Executing $(cmd)")
        r=run(cmd)
        if r.exitcode != 0
                println("Exit Code: $(r.exitcode)")
                println("Error: $(r.err)")
                @error("make failed $(cmd)")
        end

        # install
        cmd=`make install`
        @info("Executing $(cmd)")
        r=run(cmd)
        if r.exitcode != 0
                println("Exit Code: $(r.exitcode)")
                println("Error: $(r.err)")
                @error("make failed $(cmd)")
        end

end

function configure_pip(o::Arguments)
        cd("$(o.destination)")

        cmd=`bin/pip3 install --upgrade pip`
        @info("Executing $(cmd)")
        r=run(cmd)
        if r.exitcode != 0
                println("Exit Code: $(r.exitcode)")
                println("Error: $(r.err)")
                @error("make failed $(cmd)")
        end

        cmd=`bin/pip3 install wheel`
        @info("Executing $(cmd)")
        r=run(cmd)
        if r.exitcode != 0
                println("Exit Code: $(r.exitcode)")
                println("Error: $(r.err)")
                @error("make failed $(cmd)")
        end

end

function install_packages(o::Arguments)
        cd("$(o.destination)")

        cmd=`bin/pip3 install swat saspy`
        @info("Executing $(cmd)")
        r=run(cmd)
        if r.exitcode != 0
                println("Exit Code: $(r.exitcode)")
                println("Error: $(r.err)")
                @error("pip3 failed $(cmd)")
        end

        if o.requirements !== nothing
                cmd=`bin/pip3 install -r $(o.requirements)`
                @info("Executing $(cmd)")
                r=run(cmd)
                if r.exitcode != 0
                        println("Exit Code: $(r.exitcode)")
                        println("Error: $(r.err)")
                        @error("pip3 failed $(cmd)")
                end
        end

end

function main()
        r = parsecmd()

        download_python(r)
        extract_python(r)
        build_python(r)
        configure_pip(r)
        install_packages(r)

end

main()

max_stars_repo_name: asbisen/build-python

max_issues_repo_name: asbisen/build-python

max_forks_repo_name: asbisen/build-python


text:
include("../linesearch/Linesearch.jl")

# Quasi-Newton-Method (BFGS), or Gradient-Descent, depending on parameter
function qn(x0, f, g,eps, maxiter, maxLSiter, maxzoomiter, c1, c2,ls, dir)
    @printf("iter \t\t fval \t\t\t diff \t\t opt \t\t alp \t\t lsiter \n")
    x_k = x0
    x_old = x_k
    oval = Inf
    H_k = Inf
    grad_prev = 0
    grad_new = 0
    val_array = zeros(0)
    stop_array = zeros(0)

    grad_new = g(x_k)
    for i in 1:maxiter
        fval = f(x_k)
        stoppingCriteria = norm(grad_new,Inf)

        append!( val_array, fval )
        append!( stop_array, stoppingCriteria )

        if stoppingCriteria < eps
           return (x_k,"optimal",val_array,stop_array)
        end
        if abs(oval - fval) < eps
            #return (x_k, "NO DECREASE",val_array,stop_array)
        end

        α = 1
        pk = -grad_new
        iter = 1

        # set search-direction
        if dir == "gd"
            pk = -grad_new
        elseif dir == "bfgs"
            if i == 1
                # n * n identity matrix
                H_0 = eye(size(x0,1))
                H_k = H_0
            else
              # approximate inverse hessian
              s_k = x_k - x_old
              # y_k is gradient difference
              y_k = grad_new - grad_prev
              rho = 1 / dot(s_k,y_k)
              Z = eye(size(x0,1)) - rho * y_k * s_k'
              H_k = Z'*H_k*Z  + rho*s_k *s_k'
              pk = -H_k * grad_new
            end
        end
        # compute stepsize by linesearch
        lsZiter = 0
        if ls == "bt"
            α,iter = backTrackingLS(f,g,maxLSiter,c1,x_k,grad_new,fval)
        elseif ls == "wolfe"
            α,iter,lsZiter = wolfeLineSearch(f,g,pk,maxLSiter,maxzoomiter,c1,c2,x_k,grad_new,fval)
        end
        @printf("%i \t\t %f \t\t %f \t %e \t %e \t\t %i %i \n",i,fval,(oval-fval),stoppingCriteria,α,iter,lsZiter)
        # update current point

        x_old = x_k
        x_k = x_k + α * pk

        # update values
        grad_prev = grad_new
        grad_new = g(x_k)
        oval = fval
    end
    return (x_k,"maxiter",val_array,stop_array)
end

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
#nprocs() == 1 && addprocs()

using GZip
using MNIST
using Knet
@everywhere using Compat

@everywhere logis(x) = 1 / (1 + exp(-x))
@everywhere dlogis(x) = logis(x) * (1 - logis(x))

#@everywhere logis(x) = max(0,x)
#@everywhere dlogis(x) = x>0?1:0
@everywhere rmse(x) = sqrt(mean(x.^2))

@everywhere function initParams(x,y,n_hidden,max_iter)
    global w1 = randn(n_hidden,size(x,1))
    global w2 = randn(n_hidden,n_hidden)
    global w3 = randn(size(y,1),n_hidden)
    B1 = rand(n_hidden,size(y,1))
    B2 = rand(n_hidden,size(y,1))
    errors = zeros(max_iter)
    outputs = zeros(size(y,2),max_iter)
    return (w1,w2,w3,errors,outputs,B1,B2)
end

@everywhere function ffward(x,y,w1,w2,w3)
    a1 = w1 * x
    h1 = logis.(a1)
    a2 = w2 * h1
    h2 = logis.(a2)
    ay = w3 * h2
    yhat = logis.(ay)
    err = yhat - y
    return(a1,h1,a2,h2,ay,yhat,err)
end

@everywhere function updateW!(w1,w2,w3,x,h1,h2,δa1,δa2,err;lr=0.1)
    w1 .-= δa1 * x' .* lr
    w2 .-= δa2 * h1' .* lr
    w3 .-= err * h2' .* lr
end

@everywhere function bp(dta,n_hidden,max_iter)
    (x,y)=dta[1]
    (w1,w2,w3,errors,outputs,_,_)=initParams(x,y,n_hidden,max_iter)
    for i = 1:max_iter
        toterr=0
        for (x,y) in dta
            (a1,h1,a2,h2,ay,yhat,err) = ffward(x,y,w1,w2,w3)
            δa2 = (w3' * err) .* dlogis.(a2)
            δa1 = (w2' * δa2) .* dlogis.(a1)
            updateW!(w1,w2,w3,x,h1,h2,δa1,δa2,err)
            toterr += rmse(err)
        end
        println(toterr)
        errors[i] = toterr/length(dta)
    end
    return errors
end

@everywhere function dfa(dta,n_hidden,max_iter)
    (x,y)=dta[1]
    (w1,w2,w3,errors,outputs,B1,B2)=initParams(x,y,n_hidden,max_iter)
    for i = 1:max_iter
        toterr=0
        for (x,y) in dta
            (a1,h1,a2,h2,ay,yhat,err) = ffward(x,y,w1,w2,w3)
            δa2 = (B2*err) .* dlogis.(a2)
            δa1 = (B1*err) .* dlogis.(a1)
            updateW!(w1,w2,w3,x,h1,h2,δa1,δa2,err)
            toterr += rmse(err)
        end
        println(toterr)
        errors[i] = toterr/length(dta)
    end
    return errors
end

@everywhere function fa(x,y,n_hidden,max_iter)
    (w1,w2,w3,errors,outputs,B1,B2)=initParams(x,y,n_hidden,max_iter)
    for i = 1:max_iter
        (a1,h1,a2,h2,ay,yhat,err) = ffward(x,y,w1,w2,w3)
        δa2 = (B2.*err) .* dlogis(a2)
        δa1 = (B1.*δa2) .* dlogis(a1)
        updateW!(w1,w2,w3,x,h1,h2,δa1,δa2,err)
        errors[i] = sum(abs(err))
        outputs[:,i] = yhat
    end
    return (outputs,errors)
end

@everywhere function ifa(x,y,n_hidden,max_iter)
    (w1,w2,w3,errors,outputs,B1,_)=initParams(x,y,n_hidden,max_iter)
    for i = 1:max_iter
        (a1,h1,a2,h2,ay,yhat,err) = ffward(x,y,w1,w2,w3)
        δa1 = (B1.*err) .* dlogis(a1)
        δa2 = (w2 * δa1) .* dlogis(a2)
        updateW!(w1,w2,w3,x,h1,h2,δa1,δa2,err)
        errors[i] = sum(abs(err))
        outputs[:,i] = yhat
    end
    return (outputs,errors)
end

@everywhere function loaddata()
    global xtrn,ytrn,xtst,ytst
    info("Loading MNIST...")
    xtrn = gzload("train-images-idx3-ubyte.gz")[17:end]
    xtst = gzload("t10k-images-idx3-ubyte.gz")[17:end]
    ytrn = gzload("train-labels-idx1-ubyte.gz")[9:end]
    ytst = gzload("t10k-labels-idx1-ubyte.gz")[9:end]
end

@everywhere function gzload(file; path=Knet.dir("data",file), url="http://yann.lecun.com/exdb/mnist/$file")
    isfile(path) || download(url, path)
    f = gzopen(path)
    a = @compat read(f)
    close(f)
    return(a)
end

@everywhere function minibatch(x, y, batchsize; atype=Array{Float32}, xrows=784, yrows=10, xscale=255)
    xbatch(a)=convert(atype, reshape(a./xscale, xrows, div(length(a),xrows)))
    ybatch(a)=(a[a.==0]=10; convert(atype, sparse(convert(Vector{Int},a),1:length(a),one(eltype(a)),yrows,length(a))))
    xcols = div(length(x),xrows)
    xcols == length(y) || throw(DimensionMismatch())
    data = Any[]
    for i=1:batchsize:xcols-batchsize+1
        j=i+batchsize-1
        push!(data, (xbatch(x[1+(i-1)*xrows:j*xrows]), ybatch(y[i:j])))
    end
    return data
end

loaddata()
global dtrn = minibatch(xtrn, ytrn, 100; atype=Array{Float32})
global dtst = minibatch(xtst, ytst, 100; atype=Array{Float32})

#x = [1 1 ; 0 1 ; 1 0 ; 0 0]'
#y = [0 1 1 0]

n_hidden = 50; # Number of hidden units
max_iter = 10; #Number of learning steps
trials = 1;

#dfa_err=@parallel (hcat) for _=1:trials
#    dfa_errors = dfa(dtrn,n_hidden,max_iter)
#    dfa_errors
#end

bp_err = @parallel (hcat) for _=1:trials
    bp_errors = bp(dtrn,n_hidden,max_iter)
    bp_errors
end

#(bp_yhat,_) = bp(x,y,n_hidden,max_iter)
#
#ifa_err=@parallel (hcat) for _=1:trials
#    (_,ifa_errors) = ifa(x,y,n_hidden,max_iter)
#    ifa_errors
#end
#(ifa_yhat,_) = ifa(x,y,n_hidden,max_iter);
#fa_err=@parallel (hcat) for _=1:trials
#    (_,fa_errors) = fa(x,y,n_hidden,max_iter)
#    fa_errors
#end
#(fa_yhat,_) = fa(x,y,n_hidden,max_iter);
#using Plots
#l = @layout [
#    a b
#    c d
#    e f
#]

#unicodeplots()
##init_notebook(true)
#p1=plot(dfa_yhat',label=["1,1" "0,1" "1,0" "0,0"],title="DFA")
#p2=plot(ifa_yhat',label=["1,1" "0,1" "1,0" "0,0"],title="IFA")
#p3=plot(fa_yhat',label=["1,1" "0,1" "1,0" "0,0"],title="FA")
#p4=plot(bp_yhat',label=["1,1" "0,1" "1,0" "0,0"],title="BP")
#p5=plot([dfa_err[:,1] ifa_err[:,1] fa_err[:,1] bp_err[:,1]],label=["DFA" "IFA" "FA" "BP"],title="First Trial Error")
#p6=plot([mean(dfa_err,2) mean(ifa_err,2) mean(fa_err,2) mean(bp_err,2)],label=["DFA" "IFA" "FA" "BP"],title="Mean Trials Error")
#plot(p1,p2,p3,p4,p5,p6,layout=l)
##plot(p1,p2,p3,p4,p5,p6,size=(900,900),layout=l)


text:
import LinearAlgebra: svd, mul!, diagm

"""
    LinearAlgebra.svd(N::T) where {T <: AbstractEcologicalNetwork}

SVD of a network (i.e. SVD of the adjacency matrix)
"""
LinearAlgebra.svd(N::T) where {T <: AbstractEcologicalNetwork} = svd(adjacency(N))

"""
    svd_truncated(N::BinaryNetwork, rnk::Integer=3)

Given a binary network `N` which adjacency matrix `A` is of dimension `n × m`,
`svd_truncated(A)` returns two matrices, `Left` and `Right`, with dimensions
respectively `n × rank` and `rank × m`, corresponding to the species
representation in the network `N`'s RDPG space of rank `rank`.

The singular value decomposition is computed using `LinearAlgebra`'s `svd`,
obtaining

`A = U * Diagonal(S) * V = U * Diagonal(√S) * Diagonal(√S) * V`.

The truncation preserves the first `rank` columns of `U * Diagonal(√S)` and the
first `rank` rows `Diagonal(√S) * V`.

We have that, `A ≃ Left * Right` (and the approximation is optimal in a
specified meaning).

#### References

Dalla Riva, G.V. and Stouffer, D.B., 2016. Exploring the evolutionary signature
of food webs' backbones using functional traits. Oikos, 125(4), pp.446-456.
https://doi.org/10.1111/oik.02305

"""
function rdpg(N::T, rnk::Integer=3) where {T<:BinaryNetwork}
    U, singular_values, V = svd(N)
    left_space = similar(U)
    right_space = similar(V')
    mul!(left_space, U, diagm(sqrt.(singular_values)))
    mul!(right_space, diagm(sqrt.(singular_values)), V')
    left_space = left_space[:,1:rnk]
    right_space = right_space[1:rnk,:]
    return left_space, right_space
end


text:
struct CMMatrix{T,C} <: AbstractMatrix{T}
    cc::Vector{C} # vector of conjugacy classes to fix the order
    r::Int # the index of conjugacy class
    m::Matrix{T} # cache of class coefficients

    function CMMatrix(
        cc::A,
        r::Int,
        T::Type = Int,
    ) where {C,A<:AbstractVector{C}}
        M = fill(-one(T), length(cc), length(cc))
        return new{T,C}(cc, r, M)
    end
end

Base.size(M::CMMatrix) = size(M.m)
Base.IndexStyle(::Type{<:CMMatrix}) = IndexCartesian()

function Base.getindex(M::CMMatrix, s::Integer, t::Integer)
    if isone(-M.m[s, t])
        out = one(first(first(M.cc)))

        M.m[s, :] .= 0
        r = M.r

        for g in M.cc[r]
            for h in M.cc[s]
                out = GroupsCore.mul!(out, g, h)
                for t = 1:size(M, 2)
                    if out == first(M.cc[t])
                        M.m[s, t] += 1
                        break
                    end
                end
            end
        end
    end
    return M.m[s, t]
end

conjugacy_classes(G::GroupsCore.Group) = conjugacy_classes_orbit(G)

function conjugacy_classes_orbit(G::GroupsCore.Group)
    id = one(G)
    S = gens(G)
    ordG = order(Int, G)

    cclasses = [PermutationGroups.Orbit([id], Dict(id => nothing))]
    elts_counted = 1

    for g in G
        any(ccl -> g ∈ ccl, cclasses) && continue
        ccl_g = PermutationGroups.Orbit(S, g, ^)
        elts_counted += length(ccl_g)
        push!(cclasses, ccl_g)
        elts_counted == ordG && break
    end

    @assert elts_counted == ordG
    return cclasses
end


text:
include("activate.jl")

using LinearAlgebra
using Random
using SparseArrays

n = 1000
k = 10000

s = rand(1:100)
@info "Seed $s"

# this seed produces numerical issues with Float64 with the k-sparse 100 lmo / for testing
s = 41
Random.seed!(s)

matrix = rand(n, n)
hessian = transpose(matrix) * matrix
linear = rand(n)
f(x) = dot(linear, x) + 0.5 * transpose(x) * hessian * x
function grad!(storage, x)
    return storage .= linear + hessian * x
end
L = eigmax(hessian)

#Run over the probability simplex and call LMO to get initial feasible point
lmo = FrankWolfe.ProbabilitySimplexOracle(1.0);
x00 = FrankWolfe.compute_extreme_point(lmo, zeros(n))

target_tolerance = 1e-5

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_accel_simplex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    hessian=hessian,
    emphasis=FrankWolfe.memory,
    L=L,
    accelerated=true,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_simplex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    hessian=hessian,
    emphasis=FrankWolfe.memory,
    L=L,
    accelerated=false,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_convex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    L=L,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

data = [trajectoryBCG_accel_simplex, trajectoryBCG_simplex, trajectoryBCG_convex]
label = ["BCG (accel simplex)", "BCG (simplex)", "BCG (convex)"]
FrankWolfe.plot_trajectories(data, label, xscalelog=true)

matrix = rand(n, n)
hessian = transpose(matrix) * matrix
linear = rand(n)
f(x) = dot(linear, x) + 0.5 * transpose(x) * hessian * x + 10
function grad!(storage, x)
    return storage .= linear + hessian * x
end
L = eigmax(hessian)

#Run over the K-sparse polytope
lmo = FrankWolfe.KSparseLMO(100, 100.0)
x00 = FrankWolfe.compute_extreme_point(lmo, zeros(n))

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_accel_simplex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    hessian=hessian,
    emphasis=FrankWolfe.memory,
    L=L,
    accelerated=true,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_simplex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    hessian=hessian,
    emphasis=FrankWolfe.memory,
    L=L,
    accelerated=false,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

x0 = deepcopy(x00)
x, v, primal, dual_gap, trajectoryBCG_convex = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    epsilon=target_tolerance,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    L=L,
    verbose=true,
    trajectory=true,
    K=1.00,
    weight_purge_threshold=1e-10,
)

data = [trajectoryBCG_accel_simplex, trajectoryBCG_simplex, trajectoryBCG_convex]
label = ["BCG (accel simplex)", "BCG (simplex)", "BCG (convex)"]
FrankWolfe.plot_trajectories(data, label, xscalelog=true)


text:
function test_uniform_grid_type(T, N, L)
  grid = UniformGrid(T, N, L)
  T == typeof(grid.L) && T == eltype(grid.zc) && T == eltype(grid.zf)
end

function test_uniform_grid_spacing(T, N, L)
  Δ = T(L/N)
  grid = UniformGrid(T, N, L)
  Δ == grid.Δc && Δ == grid.Δf
end

function test_uniform_grid_limits_zc_left(T, N, L)
  Δ = L/N
  grid = UniformGrid(T, N, L)
  grid.zc[1] ≈ -L+0.5Δ
end

function test_uniform_grid_limits_zc_right(T, N, L)
  Δ = L/N
  grid = UniformGrid(T, N, L)
  grid.zc[N] ≈ -0.5Δ
end

function test_uniform_grid_limits_zf_left(T, N, L)
  Δ = L/N
  grid = UniformGrid(T, N, L)
  grid.zf[1] ≈ -L
end

function test_uniform_grid_limits_zf_right(T, N, L)
  Δ = L/N
  grid = UniformGrid(T, N, L)
  grid.zf[N+1] ≈ 0
end

@testset "Grids" begin
    nz, Lz = 3, 4.2
    for T in (Float64, Float32, Float16)
        @test test_uniform_grid_type(T, nz, Lz)
        @test test_uniform_grid_spacing(T, nz, Lz)
        @test test_uniform_grid_limits_zc_right(T, nz, Lz)
        @test test_uniform_grid_limits_zf_right(T, nz, Lz)
        @test test_uniform_grid_limits_zc_left(T, nz, Lz)
        @test test_uniform_grid_limits_zf_left(T, nz, Lz)
    end
end


text:
const C_HEAD = Crayon(bold=true, foreground=:blue)
const C_DIFF = YELLOW_FG

function ShowRetest()
    println(Crayon(foreground=:yellow, negative=true)("Rerunning code..."))
end

function ShowHeader(S, (func,args,kwds))
    display(CLEAR())
    # println("ZCLEARZ")
    # run(`clear`)
    C = NEGATIVE

    println(C("Command:"), " $func(", join(args, ", "), " ; ", join(kwds, ", "), ")")
    time = Dates.format(now(), dateformat"HH:MM:SS")
    println(C("Iteration: $(lpad(S.n,4))"), " - time: $(time)")
    if S.repeats >= 2
        println(Crayon(foreground=:green)("Repeated: $(lpad(S.repeats,4))"))
    else
        println()
    end
end

function ShowError(obs::OBSERVATION)
    @assert obs.result isa Exception
    bts = Base.process_backtrace(obs.stacktrace)
    rows = map(bts) do bt
        frame = bt[1]
        fileinfo = basename(string(frame.file)) * ":$(frame.line)" * (frame.inlined ? " [inlined]" : "")
        funcinfo = sprint(Base.StackTraces.show_spec_linfo, frame)

        out = [fileinfo, funcinfo]
    end

    if isempty(rows)
        rows = []
    else
        tab = hcat(rows...) |> permutedims

        # TODO: Turn this into a generic table print. Or find where others do this.
        total_max = 80
        max_sizes = [40 40]
        pad_funcs = [lpad rpad]
        col_sizes = maximum(length.(tab), dims=1)
        col_sizes = min.(col_sizes, max_sizes)
        if sum(col_sizes) < total_max
            col_sizes[end] += total_max - sum(col_sizes)
        end

        Truncate(s,n,pad) = length(s) > n ? s[1:nextind(s,n-1)] : pad(s,n)
        tab = Truncate.(tab, col_sizes, pad_funcs)

        rows = join.(eachrow(tab), " | ")
    end
    println(RED_FG("Exception: ", sprint(showerror, obs.result)))
    println(YELLOW_FG(join(rows,"\n")))
    println()
    ShowStds(obs)
    println()
end

ShowObservation(obs::Nothing, last_success) = println("No observation!")
function ShowObservation(obs, last_success)
    # println()
    print(C_HEAD("Duration: "), PrettyTime(obs.time))
    if last_success !== nothing
        Δt = obs.time - last_success.time
        if Δt > 0
            print(RED_FG(" + ", PrettyTime(abs(Δt))))
        else
            print(GREEN_FG(" - ", PrettyTime(abs(Δt))))
        end
    end
    println()
    if (obs.result isa Exception)
        PrettyResult(last_success, nothing)
    else
        PrettyResult(obs, last_success)
    end
end

ShowStds(obs::Nothing) = nothing
function ShowStds(obs)
    if !isempty(obs.stdout)
        println()
        println(GREEN_FG("With stdout output: "))
        println(obs.stdout)
    end
    if !isempty(obs.stderr)
        println()
        println(RED_FG("With stderr output: "))
        println(obs.stderr)
    end
end

PrettyResult(obs::Nothing, last_success) = println("No observation")
function PrettyResult(obs, last_success)
    print(C_HEAD("Return: "))

    if last_success === nothing || typeof(last_success.result) == typeof(obs.result)
        type_crayon = Crayon(faint=true)
    else
        type_crayon = C_DIFF
    end
    print(type_crayon(string(typeof(obs.result))))

    if obs.inferred_type != typeof(obs.result)
        print(RED_FG(" <: " * string(obs.inferred_type)))
    end

    println()

    try
        display(obs.result)
    catch exc
        println(RED_FG("Got error while trying to display result."))
        showerror(stderr, exc)
    end
    println()

    ShowStds(obs)
end

# This is a setup to allow it to be overridden in other displays.
struct CLEAR end

import Base: show, display
show(obj::CLEAR) = println(string(obj))

display(d::AbstractDisplay, ::CLEAR) = run(`clear`)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineScaleSetProperties <: SwaggerModel
    upgradePolicy::Any # spec type: Union{ Nothing, UpgradePolicy } # spec name: upgradePolicy
    automaticRepairsPolicy::Any # spec type: Union{ Nothing, AutomaticRepairsPolicy } # spec name: automaticRepairsPolicy
    virtualMachineProfile::Any # spec type: Union{ Nothing, VirtualMachineScaleSetVMProfile } # spec name: virtualMachineProfile
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    overprovision::Any # spec type: Union{ Nothing, Bool } # spec name: overprovision
    doNotRunExtensionsOnOverprovisionedVMs::Any # spec type: Union{ Nothing, Bool } # spec name: doNotRunExtensionsOnOverprovisionedVMs
    uniqueId::Any # spec type: Union{ Nothing, String } # spec name: uniqueId
    singlePlacementGroup::Any # spec type: Union{ Nothing, Bool } # spec name: singlePlacementGroup
    zoneBalance::Any # spec type: Union{ Nothing, Bool } # spec name: zoneBalance
    platformFaultDomainCount::Any # spec type: Union{ Nothing, Int32 } # spec name: platformFaultDomainCount
    proximityPlacementGroup::Any # spec type: Union{ Nothing, SubResource } # spec name: proximityPlacementGroup
    hostGroup::Any # spec type: Union{ Nothing, SubResource } # spec name: hostGroup
    additionalCapabilities::Any # spec type: Union{ Nothing, AdditionalCapabilities } # spec name: additionalCapabilities
    scaleInPolicy::Any # spec type: Union{ Nothing, ScaleInPolicy } # spec name: scaleInPolicy

    function VirtualMachineScaleSetProperties(;upgradePolicy=nothing, automaticRepairsPolicy=nothing, virtualMachineProfile=nothing, provisioningState=nothing, overprovision=nothing, doNotRunExtensionsOnOverprovisionedVMs=nothing, uniqueId=nothing, singlePlacementGroup=nothing, zoneBalance=nothing, platformFaultDomainCount=nothing, proximityPlacementGroup=nothing, hostGroup=nothing, additionalCapabilities=nothing, scaleInPolicy=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetProperties, Symbol("upgradePolicy"), upgradePolicy)
        setfield!(o, Symbol("upgradePolicy"), upgradePolicy)
        validate_property(VirtualMachineScaleSetProperties, Symbol("automaticRepairsPolicy"), automaticRepairsPolicy)
        setfield!(o, Symbol("automaticRepairsPolicy"), automaticRepairsPolicy)
        validate_property(VirtualMachineScaleSetProperties, Symbol("virtualMachineProfile"), virtualMachineProfile)
        setfield!(o, Symbol("virtualMachineProfile"), virtualMachineProfile)
        validate_property(VirtualMachineScaleSetProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(VirtualMachineScaleSetProperties, Symbol("overprovision"), overprovision)
        setfield!(o, Symbol("overprovision"), overprovision)
        validate_property(VirtualMachineScaleSetProperties, Symbol("doNotRunExtensionsOnOverprovisionedVMs"), doNotRunExtensionsOnOverprovisionedVMs)
        setfield!(o, Symbol("doNotRunExtensionsOnOverprovisionedVMs"), doNotRunExtensionsOnOverprovisionedVMs)
        validate_property(VirtualMachineScaleSetProperties, Symbol("uniqueId"), uniqueId)
        setfield!(o, Symbol("uniqueId"), uniqueId)
        validate_property(VirtualMachineScaleSetProperties, Symbol("singlePlacementGroup"), singlePlacementGroup)
        setfield!(o, Symbol("singlePlacementGroup"), singlePlacementGroup)
        validate_property(VirtualMachineScaleSetProperties, Symbol("zoneBalance"), zoneBalance)
        setfield!(o, Symbol("zoneBalance"), zoneBalance)
        validate_property(VirtualMachineScaleSetProperties, Symbol("platformFaultDomainCount"), platformFaultDomainCount)
        setfield!(o, Symbol("platformFaultDomainCount"), platformFaultDomainCount)
        validate_property(VirtualMachineScaleSetProperties, Symbol("proximityPlacementGroup"), proximityPlacementGroup)
        setfield!(o, Symbol("proximityPlacementGroup"), proximityPlacementGroup)
        validate_property(VirtualMachineScaleSetProperties, Symbol("hostGroup"), hostGroup)
        setfield!(o, Symbol("hostGroup"), hostGroup)
        validate_property(VirtualMachineScaleSetProperties, Symbol("additionalCapabilities"), additionalCapabilities)
        setfield!(o, Symbol("additionalCapabilities"), additionalCapabilities)
        validate_property(VirtualMachineScaleSetProperties, Symbol("scaleInPolicy"), scaleInPolicy)
        setfield!(o, Symbol("scaleInPolicy"), scaleInPolicy)
        o
    end
end # type VirtualMachineScaleSetProperties

const _property_map_VirtualMachineScaleSetProperties = Dict{Symbol,Symbol}(Symbol("upgradePolicy")=>Symbol("upgradePolicy"), Symbol("automaticRepairsPolicy")=>Symbol("automaticRepairsPolicy"), Symbol("virtualMachineProfile")=>Symbol("virtualMachineProfile"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("overprovision")=>Symbol("overprovision"), Symbol("doNotRunExtensionsOnOverprovisionedVMs")=>Symbol("doNotRunExtensionsOnOverprovisionedVMs"), Symbol("uniqueId")=>Symbol("uniqueId"), Symbol("singlePlacementGroup")=>Symbol("singlePlacementGroup"), Symbol("zoneBalance")=>Symbol("zoneBalance"), Symbol("platformFaultDomainCount")=>Symbol("platformFaultDomainCount"), Symbol("proximityPlacementGroup")=>Symbol("proximityPlacementGroup"), Symbol("hostGroup")=>Symbol("hostGroup"), Symbol("additionalCapabilities")=>Symbol("additionalCapabilities"), Symbol("scaleInPolicy")=>Symbol("scaleInPolicy"))
const _property_types_VirtualMachineScaleSetProperties = Dict{Symbol,String}(Symbol("upgradePolicy")=>"UpgradePolicy", Symbol("automaticRepairsPolicy")=>"AutomaticRepairsPolicy", Symbol("virtualMachineProfile")=>"VirtualMachineScaleSetVMProfile", Symbol("provisioningState")=>"String", Symbol("overprovision")=>"Bool", Symbol("doNotRunExtensionsOnOverprovisionedVMs")=>"Bool", Symbol("uniqueId")=>"String", Symbol("singlePlacementGroup")=>"Bool", Symbol("zoneBalance")=>"Bool", Symbol("platformFaultDomainCount")=>"Int32", Symbol("proximityPlacementGroup")=>"SubResource", Symbol("hostGroup")=>"SubResource", Symbol("additionalCapabilities")=>"AdditionalCapabilities", Symbol("scaleInPolicy")=>"ScaleInPolicy")
Base.propertynames(::Type{ VirtualMachineScaleSetProperties }) = collect(keys(_property_map_VirtualMachineScaleSetProperties))
Swagger.property_type(::Type{ VirtualMachineScaleSetProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetProperties }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetProperties[property_name]

function check_required(o::VirtualMachineScaleSetProperties)
    true
end

function validate_property(::Type{ VirtualMachineScaleSetProperties }, name::Symbol, val)
end


text:
using CairoMakie, LinearAlgebra
CairoMakie.activate!()
function E(q, rx, ry, x, y)
    d = sqrt((x-rx)^2 + (y-ry)^2)^3
    return (q * (x - rx) / d, q * (y - ry) / d)
end
function charges(; nq = 2)
    qs = []
    for i in 1:nq
        q = i % 2 * 2 - 1
        push!(qs, (q, cos(2π * i / nq), sin(2π * i / nq)))
    end
    qs
end
function fieldE(x,y)
    Ex, Ey = 0, 0
    for q in qs
        ex, ey = E(q..., x, y)
        Ex += ex
        Ey += ey
    end
    Point(Ex, Ey)
end
set_theme!(theme_dark())
fig = Figure(resolution = (600,400))
ax1 = Axis(fig[1,1]; aspect = DataAspect())
ax2 = Axis(fig[1,2]; aspect = DataAspect())
# bug alert, arrows are still not quite in the right direction
qs = charges()
streamplot!(ax1, fieldE, -2..2, -2..2; arrow_size = 6, linewidth = 0.5,
    colorrange = (-3,3), colormap = :CMRmap)
[scatter!(ax1, Point(qs[i][2:3]), color = qs[i][1] > 0 ? :red : :dodgerblue,
    markersize = 10) for i in 1:length(qs)]
qs = charges(; nq = 4)
streamplot!(ax2, fieldE, -2..2, -2..2; arrow_size = 6, linewidth = 0.5,
    colorrange = (-3,3), colormap = :CMRmap)
[scatter!(ax2, Point(qs[i][2:3]), color = qs[i][1] > 0 ? :red : :dodgerblue,
    markersize = 10) for i in 1:length(qs)]
limits!(ax, -2,2,-2,2)
hidedecorations!(ax1; grid = false)
hidedecorations!(ax2; grid = false)
save("streamplot.png", fig, px_per_unit = 2)
fig


text:
__precompile__(true)

module AutoViz

using Reexport
using Parameters
using AutomotiveDrivingModels

@reexport using Colors
@reexport using Cairo

import Reel
Reel.set_output_type("gif")

# using Reactive
# import Gtk
# using PGFPlots

export
        DEFAULT_CANVAS_WIDTH,
        DEFAULT_CANVAS_HEIGHT,
        render!,
        get_pastel_car_colors

const DEFAULT_CANVAS_WIDTH = 1000
const DEFAULT_CANVAS_HEIGHT = 600

include("colorscheme.jl")
include("rendermodels.jl")

include("cameras.jl")
include("interface.jl")
include("overlays.jl")
include("reel_drive.jl")

include("1d/main.jl")
include("2d/main.jl")

end # module


text:
using UnPack: @unpack

using .FiniteStrains: ToEulerianStrain, ToNaturalStrain, EulerianStrain, NaturalStrain

import .FiniteStrains: straintype

export getparam, orderof

include("types.jl")
include("ev.jl")
include("pv.jl")
include("bv.jl")

_ispositive(x) = x > zero(x)  # Do not export!

# Ref: https://github.com/JuliaLang/julia/blob/4a2830a/base/array.jl#L125
"""
    orderof(x::FiniteStrainParameters)

Return the order of a `FiniteStrainParameters`.

# Examples
```jldoctest
julia> orderof(BirchMurnaghan(40, 0.5, 4, 0)) == 3
true
```
"""
orderof(::Type{<:FiniteStrainParameters{N}}) where {N} = N
orderof(x::FiniteStrainParameters) = orderof(typeof(x))

function Base.show(io::IO, param::Parameters)  # Ref: https://github.com/mauro3/Parameters.jl/blob/3c1d72b/src/Parameters.jl#L542-L549
    if get(io, :compact, false)
        Base.show_default(IOContext(io, :limit => true), param)
    else
        # just dumping seems to give ok output, in particular for big data-sets:
        T = typeof(param)
        println(io, T)
        for f in propertynames(param)
            println(io, " ", f, " = ", getproperty(param, f))
        end
    end
end

"""
    getparam(eos::EquationOfStateOfSolids)

Get the `Parameters` from an `EquationOfStateOfSolids`.
"""
getparam(eos::EquationOfStateOfSolids) = eos.param

straintype(::Type{<:BirchMurnaghan}) = EulerianStrain
straintype(::Type{<:PoirierTarantola}) = NaturalStrain
straintype(x::FiniteStrainParameters) = straintype(typeof(x))

Base.eltype(::Type{<:Parameters{T}}) where {T} = T


text:
using AWS.EC2
using AWS

include("config.jl")

n=2
## env=AWSEnv(; id=id, key=key, dbg=dbg, region=region)
env=AWSEnv()
instances = ec2_launch(ami, keyname; env=env, owner=owner, insttype=insttype, n=n, clustername="ec2test")
println(instances)

all_instances=ec2_instances_by_tag("Owner", owner, env=env, running_only=false)
println(all_instances)

try
    ec2_show_status(instances)

    hostnames = ec2_hostnames(instances)
    println(hostnames)

    newp = ec2_addprocs(instances, keyfile; dir=dir_on_ami, exename=joinpath(dir_on_ami, "julia"))
    if length(newp) != (n*2)
        println("ERROR : Did not auto-detect number of cores $(length(newp)) != $(2*n)")
    end
#    println("Launched workers on instance type $insttype pids: \n$newp ")

    newp = ec2_addprocs(instances, keyfile; dir=dir_on_ami, exename=joinpath(dir_on_ami, "julia"), workers_per_instance=3)
    if length(newp) != (n*3)
        println("ERROR : Did not launch required number of workers $(length(newp)) != $(3*n)")
    end
#    println("Launched workers with pids: \n$newp")

    newp = ec2_addprocs(instances, keyfile; dir=dir_on_ami, exename=joinpath(dir_on_ami, "julia"), num_workers=5)
    if length(newp) != 5
        println("ERROR : Did not launch required number of workers $(length(newp)) != 5")
    end
#    println("Launched workers with pids: \n$newp")

    newp = ec2_addprocs(instances, keyfile; dir=dir_on_ami, exename=joinpath(dir_on_ami, "julia"), num_workers=1)
    if length(newp) != 1
        println("ERROR : Did not launch required number of workers $(length(newp)) != 1")
    end
#    println("Launched workers with pids: \n$newp")

catch e
    println("ERROR: $e")
    # We anyway need to terminate the instances....
end

resp = EC2.DescribeInstances(env)
## @show resp

EC2.check_running(env, instances)
EC2.MonitorInstances(env)

rmprocs(workers())

println("Sleeping for 20 secs !!!")
sleep(20)

println("Stopping instances !!!")
ec2_stop(instances; env=env)

println("Terminating instances !!!")
@show instances
ec2_terminate(instances; env=env)

println("Checking instance status !!!")
ec2_show_status(instances; env=env)


text:
using Test
using ProjectEuler

@testset "ProjectEuler" begin
    @test pe1()==233168
    @test pe2()==4613732
    @test pe3()==6857
    @test pe4()==906609
    @test pe5()==232792560
    @test pe6()==25164150
    @test pe7()==104743
    @test pe8()==23514624000
    @test pe9()==31875000
    @test pe10()==142913828922
    @test pe11()==70600674
    @test pe12()==76576500
    @test pe13()==5537376230
    @test pe14()==837799
    @test pe15()==137846528820
    @test pe16()==1366
    @test pe17()==21124
    @test pe18()==1074

    @test pe20()==648

    @test pe22()==871198282

    @test pe25()==4782

    @test pe34()==40730
end

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_name: miRoox/ProjectEuler