text:
function isEven1(a = 1, b = 10)
    try
        if a >= 0 && b >= 0
            distance = sumInterval(a, b) #1. distance
            _isPositive = isPositive(distance) #2. positive distance
            if _isPositive #positive?
                isEven = distance % 2 == 0 ? true : false
                if isEven           # even?
                    return true
                elseif !isEven # || isPositive
                    return false
                else
                    println(UnexpMsg)
                end
            else
                throw(error("Positive number error"))
            end
            #end
        end
    catch positiveError
        #Exception
        @error "ERROR: " * positiveMsg exception = (positiveError, catch_backtrace())
        println(positiveMsg) #positive arguments error
    end
end
#=error handling  https://discourse.julialang.org/t/how-to-find-the-error-line-using-try-catch-statements/65511/2
credits: @skoffer =#
#=1. either:
1.  you rethrow error, OR
2. you can use `catch_backtrace()` with any of logger macros

me: playing frisbe with try catch using catch_backtrace: as easy as 1,2,3  =#
try
    condition = true  #1.check condition
    if condition # if it works, proceed
    else
        throw(error("Positive number error")) # 2. throw(error(ExceptionError))
    end

catch positiveError # 3. catch (positiveError object )
    @error "ERROR: " * positiveMsg exception = (positiveError, catch_backtrace())   # define Exception here, passing arguments 1. positiveError object, 2. call catch_backtrace() (to catch it)
end

#end
#=
isEven1(-1, 10)
isEven1(1, -10)
isEven1(-1, -10)
benchmark isEven1()=#

#--- isEven2

@propagate_inbounds function isEven2(a = 1, b = 10) # = readable #preferred #fast  #optimized  #2 if-statements
    try
        @inbounds if a >= 0 && b >= 0   #==#
            sumInt = sumInterval(a, b) #1.distance
            evenCond = sumInt % 2 == 0 ? true : false
            _isPositive = isPositive(sumInt)
            @inbounds if evenCond && _isPositive  ? true : false  #ok
        end
        else
            throw(error("Positive number error"))
        end
    catch positiveError
        #Exception
        @error "ERROR: " * positiveMsg exception = (positiveError, catch_backtrace())
        println(positiveMsg) #positive arguments error
    end #==#
end

#=
isEven1(-1, 10)
isEven1(1, -10)
isEven1(-1, -10)=#

@benchmark isEven2(10^7)

@propagate_inbounds function isEven2(m) # = readable #preferred #fast  #optimized  #2 if-statements
    try
        @inbounds if m >= 0  # && b >= 0   #==#
            distance = sumInterval(m) #1.distance
            evenCond = distance % 2 == 0 ? true : false
            _isPositive = isPositive(distance)
            @inbounds if evenCond && _isPositive #ok ? true : false
            end
        else
            throw(error("Positive number error"))
        end
    catch positiveError
        #Exception
        @error "ERROR: " * positiveMsg exception = (positiveError, catch_backtrace())
        println(positiveMsg) #positive arguments error
    end #==#
end

@inline @propagate_inbounds function replaceVector(v = [2, 3], a = [1, 2, 3, 4]; i = 1)
    lenV = length(v)
    @inbounds if lenV < length(a) # first assumption
        @inbounds a[i:lenV] = v[i:lenV]

    else
        println(UnexpMsg)
    end
    return a
end
#=
BenchmarkTools.Trial: 10000 samples with 960 evaluations.
 Range (min … max):   86.979 ns …   4.819 μs  ┊ GC (min … max):  0.00% … 97.34%
 Time  (median):      95.729 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   126.652 ns ± 201.347 ns  ┊ GC (mean ± σ):  10.76% ±  6.80%

  ▇██▅▄▄▄▄▄▃▃▃▃▃▂▂▂▂▂▂▁▂▁▁▁▁▁   ▁▁                              ▂
  ███████████████████████████████████████▇▇▇█▇▇▇▅▆▅▆▁▆▅▆▅▅▃▅▅▄▅ █
  87 ns         Histogram: log(frequency) by time        302 ns <

 Memory estimate: 256 bytes, allocs estimate: 3.

 =#

@inline function replaceVector2(v = [2, 3], a = [1, 2, 3, 4]; i = 1)
    lenV = length(v)
    @inbounds if lenV < length(a) # first assumption
        @inbounds a[i:lenV] = v[i:lenV]

    else
        println(UnexpMsg)
    end
    return a
end

#=
BenchmarkTools.Trial: 10000 samples with 955 evaluations.
 Range (min … max):   87.016 ns …  13.450 μs  ┊ GC (min … max):  0.00% … 98.67%
 Time  (median):      95.183 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   132.741 ns ± 239.163 ns  ┊ GC (mean ± σ):  11.10% ±  6.89%

  ▇█▅▄▃▃▄▃▃▃▃▂▂▂▂▂▂▂▁▁▁▁▁▂▁▁▁▁                                  ▂
  █████████████████████████████████▇▇▆▇▇▇▇▇▆▆▇▆▆▅▅▅▆▆▅▆▅▅▅▄▄▄▁▅ █
  87 ns         Histogram: log(frequency) by time        369 ns <

 Memory estimate: 256 bytes, allocs estimate: 3.

###################

 =#

 @noninline function replaceVector3(v = [2, 3], a = [1, 2, 3, 4]; i = 1)
    lenV = length(v)
    @inbounds if lenV < length(a) # first assumption
        @inbounds a[i:lenV] = v[i:lenV]

    else
        println(UnexpMsg)
    end
    return a
end

@noinline @propagate_inbounds function replaceVector4(v = [2, 3], a = [1, 2, 3, 4]; i = 1)
    lenV = length(v)
    @inbounds if lenV < length(a) # first assumption
        @inbounds a[i:lenV] = v[i:lenV]

    else
        println(UnexpMsg)
    end
    return a
end
#=
BenchmarkTools.Trial: 10000 samples with 960 evaluations.
 Range (min … max):   86.979 ns …   4.819 μs  ┊ GC (min … max):  0.00% … 97.34%
 Time  (median):      95.729 ns               ┊ GC (median):     0.00%
 Time  (mean ± σ):   126.652 ns ± 201.347 ns  ┊ GC (mean ± σ):  10.76% ±  6.80%

  ▇██▅▄▄▄▄▄▃▃▃▃▃▂▂▂▂▂▂▁▂▁▁▁▁▁   ▁▁                              ▂
  ███████████████████████████████████████▇▇▇█▇▇▇▅▆▅▆▁▆▅▆▅▅▃▅▅▄▅ █
  87 ns         Histogram: log(frequency) by time        302 ns <

 Memory estimate: 256 bytes, allocs estimate: 3.
=#

#inline functions
@benchmark replaceVector  #noinline better
#=
BenchmarkTools.Trial: 10000 samples with 1000 evaluations.
 Range (min … max):  0.001 ns … 0.100 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     0.001 ns             ┊ GC (median):    0.00%
 Time  (mean ± σ):   0.035 ns ± 0.047 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  █                                                       ▁
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▂
  0.001 ns       Histogram: frequency by time        0.1 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

@benchmark replaceVector2

#no inline
@benchmark replaceVector3
@benchmark replaceVector4

#-----------------
res = compareVector() #<------          Attention!
α = 1;
β = 2;
arr = [2, 1, 3, 4]
buildRangeAroundPoint(α
, mid, β)
println(res)
typeof(res)

#-----------------
@propagate_inbounds function buildRangeAroundPoint(a, mid, b) #checked  # the point of buildingRanges is the one I'm concerned about
    if a >= 0 && mid >= 0 && b >= 0
        q = []
        @inbounds push!(q, buildInterval(a, mid))
        @inbounds push!(q, buildInterval(mid + 1, b))

        return q
    else
        println(positiveMsg)
    end
end
v = collect((2, 2))
rr = [1, 2, 3, 4]
res = replaceVector(v, rr) #done

#--later
@propagate_inbounds reduce(vcat, map(exp, reduce(hcat, transpose(res))))  #try maping data (vector) to a function of choice (exp)

v = [2, 2]
a = [1, 2, 3, 4]

#-------

#-----------------

permutedims(reduce(hcat, (collect(1:3) for i = 1:2)))#only hcat
reduce(vcat, transpose([collect(1:3) for i = 1:2])) # vcat, with transpose #ok
permutedims(reduce(vcat, transpose([collect(1:3) for i = 1:2])))
reduce(vcat, transpose([collect(1:3) for i = 1:2]))#ERROR: no vcat  # done [;,;,1] #above line-transpose
stack(collect(1:3) for i = 1:2; dims = 1) #not for julia 1.8

#--------

#--------
#---test ---
mid = 2.5
above = Int(ceil(mid))
below = Int(floor(mid))

#@propagate_inbounds
function middle(st, ed)
    try
        q = []
        cond = isEven(st, ed)
        @inbounds if cond
            #return mid , mid+1
            mid = Int(ϟ(st, ed) // 2)
            # return mid, mid + 1
            push!(q, mid)
            push!(q, mid + 1)
            #@inbounds
        elseif !cond #mid = 2.5 #inapplicable for an index
            mid = ϟ(st, ed) / 2 # floating-point division sumInterval(a, b) / 2 * 1.0 # freely allowing floats, to be ceiled & floored
            # above = Int(ceil(check)) #nearest index above
            #  below = Int(floor(check))
            above = Int(ceil(mid))
            below = Int(floor(mid))
            push!(q, below)
            push!(q, above)
            # return q  #return below, above
            # @inbounds
        else #throws erro
            throw(error(UnexpObj)) #generate error Object UnexpObj

        end # placecatch here @error $stringMsg exception = (ErrorObject, catch_backtrace())
    catch UnexpMsg
        exception = (UnexpObj, catch_backtrace)
    end   # or if anything is smooth, return as below
    return q
end# final end
#catch  the UnexpObj - on the catch_backtrace
#end
#=

    q = []
    cond = isEven(st, ed)
    @inbounds if cond
        #return mid , mid+1
        mid = Int(ϟ(st, ed) // 2)
        # return mid, mid + 1
        push!(q, mid)
        push!(q, mid + 1)
        #@inbounds
    elseif !cond #mid = 2.5 #inapplicable for an index
        mid = ϟ(st, ed) / 2 # floating-point division sumInterval(a, b) / 2 * 1.0 # freely allowing floats, to be ceiled & floored
        # above = Int(ceil(check)) #nearest index above
        #  below = Int(floor(check))
        above = Int(ceil(mid))
        below = Int(floor(mid))
        push!(q, below)
        push!(q, above)
        # return q  #return below, above
        # @inbounds
    else #throws erro  #generate error Object UnexpObj
        throw(error(UnexpObj))

    end
    return q
end# final end - placecatch here @error $stringMsg exception = (ErrorObject, catch_backtrace())
catch
@error UnexpMsg exception =(UnexpObj, catch_backtrace())
println(positiveMsg) #positive arguments error
end
=#
_midV = middle(1, 10) # yes, mid = 5 return 5,6

#--- we are here [current bottleNeck ]
#requires buildInterval
#uses triad if-elseif-else
#TODO:solution

#Don't run this codeblock!
#=
@propagate_inbounds function doCompare(st = 1, ed = 2, a = [2, 1, 3, 4]) #errrorneous # a bit absurd (don't you think?) #no everthing is fine # it's your nagging mind that is
    try
        _first = copy(a[st])
        _last = copy(a[ed])  #creates a shallow copy (what's desired for optimization)
        @inbounds if a[st] > a[ed] #valueAt(2) > valueAt(1) isa true  #_first > _last #
            @inbounds a[st], a[ed] = a[ed], a[st]        #an inbounds swap #actual array swap
        elseif _first < _last # <--------- necessary step, to jump later into errorchecking else
        else #can throw error - here
            throw()
            println(UnexpMsg)
        end
        v = buildInterval(_first, _last)  #collect((_first, _last))
        return v # ERROR: BoundsError: attempt to access 2-element Vector{Int64} at index [10]
    catch

    end

end
doCompare(1, 2)=#
#---------Start from here
#issue reconstruction
a = [2, 1, 3, 4]
st = 1;
ed = 2; # todo move st ed until end of arr

#1. check comparison  #check
#2. check _first & _last   - a copy only copies the value of content (lightweight) - Desired  #check

_first = copy(a[st])
_last = copy(a[ed])

#works fine  -fair enough # add only this condition - which we're intereste to evaluate

@inbounds if _first > _last #valueAt(2) > valueAt(1) isa true  #_first > _last #
    #Base.@propagate_inbounds# directly change contents of the Original Array
    @inbounds a[st], a[ed] = a[ed], a[st]        #an inbounds swap #actual array swap
end
#elseif  _first < _last #this line should be omitted at all costs

#skip line
#do nothing
#end

#no issues here
a #array #correct!

#---not here: apply try -catch playing safe
#should work fine, as well

try
    _first = copy(a[st])
    _last = copy(a[ed])

    #works fine  -fair enough - goal: walk me through code until we reach the catch part

    @inbounds if _first > _last #valueAt(2) > valueAt(1) isa true  #_first > _last #
        #Base.@propagate_inbounds# directly change contents of the Original Array
        println(a[st], a[ed])
        return @inbounds a[st], a[ed] = a[ed], a[st]        #an inbounds swap #actual array swap
    #end
    elseif _first < _last #possible - correct situation (to deal with)
        #Intent: skip
        return  # this  makes things work !
    else #2. throw frisbe error here
        throw(error("Unexpected Error")) # 2. throw(error(ExceptionError))
    end

catch UnexpectedError # 3. catch (UnexpectedError object )
    @error "ERROR: " * positiveMsg exception = (UnexpectedError, catch_backtrace())   # define Exception here, passing arguments 1. positiveError object, 2. call catch_backtrace() (to catch it)
end #ends try - finally afterthat return whatever correct value you've been working on  (if not already )

@propagate_inbounds function indexOf(i, v::Vector)
    try
        res = findfirst(isequal(i), v)
        typeof(res) == Nothing ? res = -1 : return Int(res) #res[1] #
    catch
        return -1
    end

end

res = findfirst(isequal(i), v)
typeof(res) == Nothing ? res = -1 : return Int(res)
#=compares vector a, it's element at first index α with second element at index β

```input:
α: first index of comparison
β: second index of comparison
a: original vector array
```

```output:
an ordered tuple of the corrected indicies (of the vector array)

```
=#

#--- compareVector
#=requires the use of indexOF,elementAt
 a = [2, 1, 3, 4]
 α = 1
 β = 2

firstContent =Int(findfirst(isequal(α), a)) #indexOf(first)
lastContent =  Int(findfirst(isequal(β), a)) #indexOf(last)

=#

#---correct
if firstContent > lastContent #correct
    a[α], a[β] = oldschoolSwap!(a[α], a[β]) #plain content swap in julia
end
a #array values are changed

function oldschoolSwap!(x, y)
    tmp = x
    x = y
    y = tmp
    return x, y
end

function compareVector(α = 1, β = 2, a = [2, 1, 3, 4])
    response = nothing
    try #1. we call this function when we'd like to compare index α with index β of a Vector array  # do your thing
        firstContent = Int(findfirst(isequal(α), a)) #indexOf(first)
        lastContent = Int(findfirst(isequal(β), a)) #indexOf(last)

        if firstContent > lastContent # correct
            response = @inbounds a[α], a[β] = oldschoolSwap!(a[α], a[β]) #plain content swap in julia  #swap array contents directly

        elseif firstContent < lastContent #only possible - correct situation (to deal with)
            #Intent: skip
            return
        else #2. throw frisbe error here
            throw(error("Unexpected Error")) # 2. throw(error(ExceptionError))
        end

    catch UnexpectedError # 3. catch `materialize` (UnexpectedError object )
        @error UnexpMsg exception = (UnexpectedError, catch_backtrace())   # define Exception here, passing arguments 1. positiveError object, 2. call catch_backtrace() (to catch it)
    end #ends try - finally afterthat return whatever correct value you've been working on  (if not already )
    return response
end

@propagate_inbounds function replaceVector(v = [2, 3], a = [1, 2, 4, 5]; i = 1)
    try
        lenV = copy(length(v))
        lenA = copy(length(a))

        @inbounds if lenV < lenA # first assumption  # |v| < |a|
            @inbounds a[i:lenV] = v[i:lenV]

        else #2. throw frisbe error here
            throw(error("Unexpected Error")) # 2. throw(error(ExceptionError))
        end

    catch UnexpectedError # 3. catch `materialize` (UnexpectedError object )
        @error UnexpMsg exception = (UnexpectedError, catch_backtrace())   # define Exception here, passing arguments 1. positiveError object, 2. call catch_backtrace() (to catch it)
    end #ends try - finally afterthat return whatever correct value you've been working on  (if not already )
    return response
end

replaceVector()

@benchmark replaceVector()

#---------------- test
v = compareVector()
typeof(compareVector)
tuple = compareVector()
v = buildInterval(tuple)# pass-in a tuple   # buildInterval(tuple[1],tuple[2])

#ends
#=
else
        println(UnexpMsg)
    end
    return a
end
=#
replaceVector() # correct
@benchmark replaceVector() # optimizable - can do better

#----

res = compareVector() #(1,2) tuple
typeof(res)#tuple
v = buildInterval(res) #to vector
a = replaceVector(v, a)
a #check a's contents

#---test indexOf
a = [2, 1, 3, 4]
indexOf(1, a)
indexOf(2, a)
@benchmark indexOf(1, a) #possible errrorneous

#=
BenchmarkTools.Trial: 96 samples with 1 evaluation.
 Range (min … max):  41.404 ms … 103.842 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     48.335 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   52.409 ms ±  12.069 ms  ┊ GC (mean ± σ):  0.49% ± 2.38%

     █▂▄█
  ▇█▆██████▇▆▆▇▆▁▅▆▆▁▆▃▁▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▁▁▁▁▃▁▁▁▁▁▁▃▃▁▁▁▁▁▁▁▃ ▁
  41.4 ms         Histogram: frequency by time          100 ms <

 Memory estimate: 1.11 MiB, allocs estimate: 14363.

indexOf (generic function with 1 method)

1

2

4-element Vector{Int64}:
 2
 1
 3
 4

2
 =#

function compareVectorExperimental([1:2;], a = [2, 1, 3, 4]) end
function compareVectorExperimental2(a = [2, 1, 3, 4], [1:2;]) end
function compareVector(α = 1, β = 2, a = [2, 1, 3, 4])

    try #1. we call this function when we'd like to compare index α with index β of a Vector array  # do your thing
        _first = Int(indexOf(α, a)) # copy(a[st])
        _last = Int(indexOf(β, a)) # copy(a[ed])
        firstIndex = a[_first]
        lastIndex = a[_last]
        @inbounds if firstIndex > lastIndex #valueAt(2) > valueAt(1) isa true  #_first > _last #
            println(a[_first], a[_last]) # debugging purposes only
            return @inbounds a[_first], a[_last] = a[_last], a[_first]     #an inbounds swap #actual array swap

        elseif firstIndex < lastIndex #only possible - correct situation (to deal with)
            #Intent: skip
            return
        else #2. throw frisbe error here
            throw(error("Unexpected Error")) # 2. throw(error(ExceptionError))
        end

    catch UnexpectedError # 3. catch `materialize` (UnexpectedError object )
        @error UnexpMsg exception = (UnexpectedError, catch_backtrace())   # define Exception here, passing arguments 1. positiveError object, 2. call catch_backtrace() (to catch it)
    end #ends try - finally afterthat return whatever correct value you've been working on  (if not already )

end

a = [2, 1, 3, 4]
_first = Int(indexOf(1, a)) # copy(a[st])
_last = Int(indexOf(2, a)) # copy(a[ed])

#------

#=
BenchmarkTools.Trial: 10000 samples with 983 evaluations.
 Range (min … max):  54.527 ns …   5.274 μs  ┊ GC (min … max): 0.00% … 97.01%
 Time  (median):     60.936 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   88.895 ns ± 126.925 ns  ┊ GC (mean ± σ):  5.20% ±  3.98%

  ▆█▅▃▂▂▂▃▂▁▃▃▃▂▃▂▂▁▁▁▁▁  ▁                                    ▁
  ████████████████████████████████████▇▆▇▇▇▇▆▆▅▆▇▇▆▅▆▅▆▆▅▅▅▅▄▅ █
  54.5 ns       Histogram: log(frequency) by time       262 ns <

=#

#------
#@benchmark replaceVector2(v,a)
@benchmark replaceVector(v, a)
#=
 Range (min … max):  21.063 ns … 130.391 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     22.568 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   26.359 ns ±  10.667 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▅█▆▃            ▁▂▂ ▁▁ ▁ ▁                                   ▁
  █████▇▇▇▇▆▅▁▃▄▄████████████▇▆▆▆▆▆▆▅▆▇▇▇▇▇▇▇▇▇▆▅▆▆▇▅▅▆▄▅▅▅▄▅▅ █
  21.1 ns       Histogram: log(frequency) by time      75.5 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.
=#

#----test-----------------------------------------------
arr = [1, 2, 3, 4]
divideConquer(arr, arr[1], length(arr) - 1)

length(arr)
#middle(st=1, ed =4)
mid = #middle(1, 4) # ambiguous function  # StackOverflowError:
    cond = isEven(mid) # ERROR: LoadError: UndefVarError: mid not defined

#--test---------------  using Findfirst: indexOF #TODO:
#include("Findfirst.jl")
#actual size

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm

function extendlines(V::Lar.Points, EV::Lar.Cells, s=1.001)
	@assert size(V,1)==2 "extendlines() is for 2D models!"
	segments = [hcat(V[:,v1], V[:,v2]) for (v1,v2) in EV]
	n = length(EV)

	for segi in segments
		# renaming of segment vertices
		x1 = segi[1,1]; x2 = segi[1,2];
		y1 = segi[2,1]; y2 = segi[2,2];
		# middle segment point
		tx = sum(x1+x2)/2
		ty = sum(y1+y2)/2
		# composite translation and scaling transformation
		segi[1,1] = s*(x1-tx)+tx; segi[1,2] = s*(x2-tx)+tx;
		segi[2,1] = s*(y1-ty)+ty; segi[2,2] = s*(y2-ty)+ty;
	end
	V = zeros(Float64, 2, 2*length(EV))
	# new vertices for expanded line segments
	for (i,segi) in enumerate(segments)
		V[:,i] = segi[:,1]
		V[:,n+i] = segi[:,2]
	end
	EV = [[k,n+k] for k=1:length(segments)]
	return V, EV
end

n=100; t=0.6
V = zeros(Float64,2,2*n)
EV = [zeros(Int64,2) for k=1:n]

for k=1:n
	v1 = rand(Float64,2)
	v2 = rand(Float64,2)
	vm = (v1+v2)/2
	transl = rand(Float64,2)
	V[:,k] = (v1-vm)*t + transl
	V[:,n+k] = (v2-vm)*t + transl
	EV[k] = [k,n+k]
end

V = Plasm.normalize(V)
model = (V,EV)
Sigma = Lar.spaceindex(model)

Plasm.view(Plasm.numbering(.1)((V,[[[k] for k=1:size(V,2)], EV])))
Z,EZ = extendlines(V,EV,1.1)
Plasm.view(Plasm.numbering(.1)((Z,[[[k] for k=1:size(Z,2)], EZ])))

Plasm.view(V,EV)
model = V,EV;
model2 = Z,EZ;
W,EW = Lar.fragmentlines(model);
Z,EZ = Lar.fragmentlines(model2);

Plasm.viewexploded(W,EW)(1.2,1.2,1.2)
Plasm.viewexploded(Z,EZ)(1.2,1.2,1.2)
W,EW = Lar.fragmentlines(model);
Z,EZ = Lar.fragmentlines(model2);
U,EVs = Lar.biconnectedComponent((W,EW::Lar.Cells));
T,ETs = Lar.biconnectedComponent((Z,EZ::Lar.Cells));
# 2-connected components (H & T)

W = convert(Lar.Points, U')
T = convert(Lar.Points, T')
cop_EV = Lar.coboundary_0(EW::Lar.Cells)
cop_EW = convert(Lar.ChainOp, cop_EV)
cop_EZ = Lar.coboundary_0(EZ::Lar.Cells)
cop_EZ = convert(Lar.ChainOp, cop_EZ)
V, copEV, copFE = Lar.Arrangement.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)
Z, copEZ, copFE2 = Lar.Arrangement.planar_arrangement(T::Lar.Points, cop_EZ::Lar.ChainOp)

triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
triangulated_faces2 = Lar.triangulate2D(Z, [copEZ, copFE2])
FVs = convert(Array{Lar.Cells}, triangulated_faces)
FZs = convert(Array{Lar.Cells}, triangulated_faces2)
W = convert(Lar.Points, V')
Z = convert(Lar.Points, Z')
Plasm.viewcolor(W::Lar.Points, FVs::Array{Lar.Cells})
Plasm.viewcolor(Z::Lar.Points, FZs::Array{Lar.Cells})

model = U,EVs;
model2 = T',ETs;
Plasm.view(Plasm.lar_exploded(model)(1.2,1.2,1.2))
Plasm.view(Plasm.lar_exploded(model2)(1.2,1.2,1.2))

ext: jl


text:
module Sparse
using SparseArrays

"""
 Allocates a block sparse matrix that stores nnz non-zero entries.
 The matrix contains only zeros after allocation.

 Input:
    rows: A vector describing the number of elements in each block row the
          block matrix has.
 columns: A vector describing the number of elements in each block column the
          block matrix has.

 Output:
       Z: Sparse matrix (all zeros)

 Example:
 >> Z = block_matrix([2, 3],[2, 3])
 >> size(Z)

 ans =

      5     5
"""
function block_matrix(rows::AbstractArray, columns::AbstractArray)
        m = sum(rows);
        n = sum(columns);

        Z = spzeros(m,n);
        return Z
end

"""
 Z = block_matrix_insert(B,rows,columns,block)
 Inserts a submatrix into a block matrix at block position (i,j)

 Input:
       B: Block matrix (see block_matrix.m for a description).
    rows: Number of elements in each block row of the block matrix.
 columns: Number of elements in each block row of the block matrix.
     i,j: Block indices denoting where the submatrix should be placed.
   block: The submatrix to insert into the block matrix.
          The size of the submatrix must match the number given in rows(i)
          and columns(j).

 Example:
 >> Z = block_matrix([2, 3],[2, 3])
 >> B = block_matrix_insert(Z,[2, 3],[2, 3],1,2,ones(2,3))
 >> B

 B =

    (1,3)        1
    (2,3)        1
    (1,4)        1
    (2,4)        1
    (1,5)        1
    (2,5)        1
>> full(B)

ans =

     0     0     1     1     1
     0     0     1     1     1
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
"""
function block_matrix_insert(B::AbstractArray,
                             rows::AbstractArray,
                             columns::AbstractArray,
                             i::Int64,
                             j::Int64,
                             block::AbstractArray)

@assert(i <= length(rows),          "Index out of bounds.")
@assert(j <= length(columns),       "Index out of bounds.")
@assert(size(block,1) == rows[i],   "Dimension mismatch")
@assert(size(block,2) == columns[j],"Dimension mismatch")

        r = vcat([1],rows)
        c = vcat([1],columns)
        row_offset = sum(r[1:i])
        col_offset = sum(c[1:j])
        B[row_offset:row_offset+r[i+1]-1,col_offset:col_offset+c[j+1]-1] = block
        return B
end

function block_matrix_2x2(a11::AbstractArray,
                          a12::AbstractArray,
                          a21::AbstractArray,
                          a22::AbstractArray)
        m1 = size(a11, 1)
        n1 = size(a11, 2)
        m2 = size(a22, 1)
        n2 = size(a22, 2)
        @assert size(a12, 1) == size(a11, 1)
        @assert size(a12, 2) == size(a22, 2)
        @assert size(a21, 1) == size(a22, 1)
        @assert size(a21, 2) == size(a11, 2)
        rows = [m1, m2]
        cols = [n1, n2]
        A = block_matrix(rows, cols)
        A = block_matrix_insert(A, rows, cols, 1, 1, a11)
        A = block_matrix_insert(A, rows, cols, 1, 2, a12)
        A = block_matrix_insert(A, rows, cols, 2, 1, a21)
        A = block_matrix_insert(A, rows, cols, 2, 2, a22)
        return A
end

function diag_block_matrix_2x2(a11::AbstractArray,
                               a22::AbstractArray)
        m1 = size(a11, 1)
        n1 = size(a11, 2)
        m2 = size(a22, 1)
        n2 = size(a22, 2)
        rows = [m1, m2]
        cols = [n1, n2]
        A = block_matrix(rows, cols)
        A = block_matrix_insert(A, rows, cols, 1, 1, a11)
        A = block_matrix_insert(A, rows, cols, 2, 2, a22)
        return A
end

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
using Dates

import InsarTimeseries: _get_pairs

split1 = convert(Array{Date,1}, [])
split2 = [Date(2017, 1, 1)]
split3 = [Date(2016, 1, 1), Date(2017, 1, 1)]

@test _get_pairs(split1) == (nothing, nothing)
@test _get_pairs(split2) == [(nothing, 2017 - 01 - 01), (2017 - 01 - 01, nothing)]
@test _get_pairs(split3) == [
    (nothing, 2016 - 01 - 01),
    (2016 - 01 - 01, 2017 - 01 - 01),
    (2016 - 01 - 01, nothing),
]


text:
using Documenter, StructuredOptimization, LinearAlgebra, DSP, AbstractFFTs, FFTW, AbstractOperators

makedocs(
  modules = [StructuredOptimization],
  format = :html,
  sitename = "StructuredOptimization",
  authors = "Niccolò Antonello and Lorenzo Stella",
  pages = Any[
  "Home"                  => "index.md",
  "Quick Tutorial Guide"  => "tutorial.md",
  "Expressions"           => "expressions.md",
  "Functions"             => "functions.md",
  "Solvers"               => "solvers.md",
  "Demos"                 => "demos.md",
  ],
)

deploydocs(
  repo   = "github.com/kul-forbes/StructuredOptimization.jl.git",
  julia  = "1.0",
  osname = "linux",
  target = "build",
  deps = nothing,
  make = nothing,
)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SLATEC_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SLATEC")
JLLWrappers.@generate_main_file("SLATEC", UUID("2967fdeb-73c1-5946-ae15-4ce8b52d4346"))
end  # module SLATEC_jll


text:
module AoC2019

include("./computer.jl")
include("./graph_algorithms.jl")

function read_numbers(filename)
  input_data = open(filename) do file
    data = readline(file)
    entries = split(data, ",")
    numbers = map(x->parse(Int64, x), entries)
    return numbers
  end
end

function read_lines(filename)
  input_data = open(filename) do file
      data = readlines(file)
      return data
  end
end

function read_strings(filename)
  input_data = open(filename) do file
      data = readlines(file)
      entries = map(x -> split(x, ","), data)
      return entries
  end
end

function read_numbers_lines(filename)
  input_data = open(filename) do file
    all_numbers = Int64[]
    data = readlines(file)
    for line in data
        entries = collect(line)
        numbers = map(x->parse(Int64, x), entries)
        append!(all_numbers, numbers)
    end
    return all_numbers
  end
end

function read_chars(filename)
  input_data = open(filename) do file
    all_numbers = Int8[]
    while !eof(file)
      c = read(file, Char)
      if c == '\n'
        continue
      end
      digit = parse(Int8, c)
      push!(all_numbers, digit)
    end
    return all_numbers
  end
end

function matlab2c(image, width, height)
  return transpose(reshape(image, width, height))
end

function visualizeArray(array, mapping)
  my_string = ""
  dims = size(array)
  local arr = map(x -> mapping[x], array)
  for i in 1:dims[1]
    for j in 1:dims[2]
      my_string *= arr[i,j]
    end
    my_string *= '\n'
  end
  return my_string
end

function pseudoNormalize(x::T)::Int32 where T<:Number
  if x > 0
    return 1
  elseif x < 0
    return -1
  else
    return 0
  end
end

function arr2text(arr)
  (x_dim, y_dim) = size(arr)
  output_str = ""
  for x in 1:x_dim
    for y in 1:y_dim
      output_str *= arr[x,y]
    end
    output_str *= '\n'
  end
  return output_str[1:end-1]
end

function text2arr(input_data)::Array{Char, 2}
  line_length = findfirst(x -> (x == '\n'), input_data) - 1
  input_data_filtered = collect(filter(x -> x != '\n', input_data))
  reshaped = reshape(input_data_filtered, line_length, :)
  return permutedims(reshaped)
end

end # module

max_stars_repo_name: buotex/aoc

max_issues_repo_name: buotex/aoc

max_forks_repo_name: buotex/aoc


text:
"""
ULMFiT - LANGUAGE MODEL [Word-by-Word]
"""

using WordTokenizers   # For accesories
using InternedStrings   # For using Interned strings
using DelimitedFiles   # For reading and writing files
using Flux  # For building models
using Flux: Tracker, crossentropy, chunk
using LinearAlgebra: norm
using BSON: @save, @load  # For saving model weights
# using CuArrays  # For GPU support

cd(@__DIR__)
include("custom_layers.jl")      # importing AWD_LSTM (ASGD Weight-Dropped LSTM)
include("utils.jl")     # including some functions
include("WikiText103_DataDeps.jl")      # including WikiText-103 Corpus

# Language Model
mutable struct LanguageModel
    vocab :: Vector
    layers :: Flux.Chain

    function LanguageModel(embedding_size::Integer=400, hidLSTMSize::Integer=1150, outLSTMSize::Integer=embedding_size;
        embedDropProb::Float64 = 0.05, wordDropProb::Float64 = 0.4, hidDropProb::Float64 = 0.5, LayerDropProb::Float64 = 0.3, FinalDropProb::Float64 = 0.4)
        vocab = intern.(string.(readdlm("vocab.csv",',', header=false)[:, 1]))
        lm = new(vocab)
        de = DroppedEmbeddings(length(vocab), embedding_size, 0.1; init = (dims...) -> init_weights(0.1, dims...))
        lm.layers = Chain(
            de,
            VarDrop(wordDropProb),
            AWD_LSTM(embedding_size, hidLSTMSize, hidDropProb; init = (dims...) -> init_weights(1/hidLSTMSize, dims...)),
            VarDrop(LayerDropProb),
            AWD_LSTM(hidLSTMSize, hidLSTMSize, hidDropProb; init = (dims...) -> init_weights(1/hidLSTMSize, dims...)),
            VarDrop(LayerDropProb),
            AWD_LSTM(hidLSTMSize, outLSTMSize, hidDropProb; init = (dims...) -> init_weights(1/hidLSTMSize, dims...)),
            VarDrop(FinalDropProb),
            x -> de(x, true),
            softmax
        )
        return lm
    end
end

Flux.@treelike LanguageModel

# Forward
function forward(lm, batch)
    batch = broadcast(x -> lm.layers[1](indices(x, lm.vocab, "_unk_")), batch)
    batch = lm.layers[2:end].(batch)
    return batch
end

# loss funciton - Loss calculation with AR and TAR regulatization
function loss(lm, gen)
    H = forward(lm, take!(gen))
    Y = broadcast(x -> gpu(Flux.onehotbatch(x, lm.vocab, "_unk_")), take!(gen))
    l = sum(crossentropy.(H, Y))
    Flux.truncate!(lm.layers)
    return l
end

# Gradient Clipping
grad_clipping(g, upper_bound) = min(g, upper_bound)

# Backward
function back!(layers, l, opt, gradient_clip::Float64)
    # Applying gradient clipping
    l = Tracker.hook(x -> grad_clipping(x, gradient_clip), l)

    # Calulating gradients and weights updation
    p = get_trainable_params(layers)
    grads = Tracker.gradient(() -> l, p)
    Tracker.update!(opt, p, grads)
    return
end

# Funciton for training Language Model
function pretrain_lm!(lm::LanguageModel; batchsize::Integer=64, bptt::Integer=70, data_loader::Channel=load_wikitext_103, gradient_clip::Float64=0.25,
    base_lr=0.004, epochs::Integer=1, checkpoint_iter::Integer=5000)

    # Initializations
    opt = ADAM(base_lr, (0.7, 0.99))    # ADAM Optimizer
    gpu!.(lm.layers)

    # Pre-Training loops
    for epoch=1:epochs
        gen = Channel(x -> generator(x, loadCorpus(); batchsize = batchsize, bptt = bptt))
        num_of_batches = take!(gen) # Number of mini-batches
        T = num_of_iters-Int(floor((num_of_iters*2)/100))   # Averaging Trigger
        set_trigger!.(T, lm.layers)  # Setting triggers for AWD_LSTM layers
        for i=1:num_of_batches

            # FORWARD PROPAGATION
            l = loss(lm, gen)

            # BACK PROPAGATION
            back!(lm.layers, l, opt, gradient_clip)
            for layer in lm.layers[1:8]
                cpu!(layer)
            end

            # ASGD Step, after Triggering
            asgd_step!.(i, lm.layers)

            reset_masks!.(lm.layers)

            println("loss: $l", " iteration number: $i")

            # Saving checkpoints
            if i == checkpoint_iter save_model!(lm) end
        end
        println("\nEpoch: $epoch")
    end
end

# To save model
function save_model!(lm::LanguageModel, filepath::String="ULMFiT-LM.bson")
    weights = Tracker.data.(params(lm))
    @save filepath weights
end

# To load model
function load_model!(lm::LanguageModel, filepath::String=joinpath(datadep"pretrained-ULMFiT", "weights.bson"))
    @load filepath weights
    Flux.loadparams!(lm, weights)
end

function load_model!()
    lm = LanguageModel()
    load_model!(lm)
    return lm
end

function load_model!(filepath::String)
    lm = LanguageModel()
    load_model!(lm, filepath)
end


text:
""" Formulas (expression trees) for computer algebra.

This module is *not* an implementation of a conventional, general-purpose CAS.
There are already many highly developed computer algebra systems. The
objectives are to:

- represent formulas as expression trees, whenever that is convenient
- display formulas in conventional mathematical notation with free variables
- allow interoperation with existing computer algebra systems
"""
module MathFormulas
export Formula, head, args, first, last, to_formula, to_wiring_diagram,
  show_latex, show_latex_formula, show_sexpr

import Base: first, last, show
import Base.Meta: show_sexpr

using AutoHashEquals
using Match

using ...Catlab
import ...Meta: strip_lines
import ...Syntax: head, args, show_latex
using ...WiringDiagrams
import ...WiringDiagrams: to_wiring_diagram

using ..AlgebraicNets
import ..AlgebraicNets: Ob, Hom, compose, id, dom, codom, otimes, munit, braid,
  mcopy, delete, mmerge, create, linear, constant, wiring, evaluate

# Data types
############

""" An expression tree for computer algebra.

We call the expression trees "formulas" to avoid confusion with Julia
expressions (`Base.Expr`) or GAT expressions (`Catlab.Syntax.GATExpr`). Usually
the operations (head symbols) are interpreted as Julia functions, e.g., `:/` is
right multiplication by the matrix pseudoinverse while `:./` is the elementwise
division.
"""
@auto_hash_equals struct Formula
  head::Symbol
  args::Vector
  Formula(head::Symbol, args...) = new(head, collect(args))
  Formula(sexp::Tuple) = new(sexp[1],
    [ x isa Tuple ? Formula(x) : x for x in sexp[2:end] ])
end
head(form::Formula) = form.head
args(form::Formula) = form.args
first(form::Formula) = first(args(form))
last(form::Formula) = last(args(form))

# Conversion
############

""" Convert Julia expression to formula.

Only a subset of the Julia syntax is supported.
"""
function to_formula(expr::Expr)
  expr_to_formula(strip_lines(expr, recurse=true))
end

function expr_to_formula(expr)
  @match expr begin
    Expr(:call, [args...]) => Formula(map(expr_to_formula, args)...)
    Expr(:->, [var::Symbol, body]) => Formula(:->, var, expr_to_formula(body))
    Expr(:&&, [args...]) => Formula(:&, map(expr_to_formula, args)...)
    Expr(:||, [args...]) => Formula(:|, map(expr_to_formula, args)...)
    Expr(:block, [arg]) => expr_to_formula(arg)
    Expr(:block, _) => error("Cannot parse multiline block as formula")
    symbol::Symbol => symbol
    literal::Number => literal
    _ => error("Cannot parse expression as formula: $expr")
  end
end

""" Convert algebraic network to formula.

Assumes that the network has a single output.
"""
function to_formula(f::AlgebraicNet.Hom, vars::Vector{Symbol})
  formulas = functor((NFormula, Formulas), f)
  @assert codom(formulas).n == 1
  substitute(formulas.terms[1], Dict(zip(formulas.vars, vars)))
end

struct NFormula
  n::Int
end
struct Formulas
  terms::Vector
  vars::Vector{Symbol}
end

""" Algebraic networks realized by formulas with free variables.

These methods should only be used with `gensym`-ed variables since they assume
that any two formulas have disjoint variables.
"""
@instance AlgebraicNetTheory(NFormula, Formulas) begin
  dom(f::Formulas) = NFormula(length(f.inputs))
  codom(f::Formulas) = NFormula(length(f.terms))

  function compose(f1::Formulas, f2::Formulas)::Formulas
    replacements = Dict(zip(f2.vars, f1.terms))
    terms = [ substitute(term, replacements) for term in f2.terms ]
    Formulas(terms, f1.vars)
  end

  function id(A::NFormula)::Formulas
    vars = gensyms(A)
    Formulas(vars, vars)
  end

  munit(::Type{NFormula}) = NFormula(0)
  otimes(A::NFormula, B::NFormula) = NFormula(A.n + B.n)

  function otimes(b1::Formulas, b2::Formulas)::Formulas
    Formulas([b1.terms; b2.terms], [b1.vars; b2.vars])
  end
  function braid(A::NFormula, B::NFormula)::Formulas
    v1, v2 = gensyms(A), gensyms(B)
    Formulas([v2; v1], [v1; v2])
  end

  mcopy(A::NFormula) = mcopy(A, 2)
  mmerge(A::NFormula) = mmerge(A, 2)
  delete(A::NFormula) = Formulas([], gensyms(A))
  create(A::NFormula) = Formulas(zeros(Int, A.n), [])

  function linear(value::Any, A::NFormula, B::NFormula)::Formulas
    nin, nout = A.n, B.n
    @assert nin == 1 && nout == 1 # FIXME: Relax this.
    var = gensym()
    Formulas([Formula(:*, value, var)], [var])
  end

  function constant(x::Any, A::NFormula)::Formulas
    @assert A.n == 1
    Formulas([x], Symbol[])
  end

  function wiring(f::Any, A::NFormula, B::NFormula)::Formulas
    vars = gensyms(A)
    terms = [ [] for i in 1:B.n ]
    for (src, tgts) in (f::WiringLayer).wires
      x = vars[src]
      for (tgt, c) in tgts
        push!(terms[tgt], c == 1 ? x : Formula(:*, c, x))
      end
    end
    Formulas(map(sum_formulas, terms), vars)
  end
end

function mcopy(A::NFormula, n::Int)::Formulas
  vars = gensyms(A)
  Formulas(reduce(vcat, fill(vars, n)), vars)
end

function mmerge(A::NFormula, n::Int)::Formulas
  vars = gensyms(A.n * n)
  Formulas([sum_formulas(vars[i:A.n:end]) for i in 1:A.n], vars)
end

Ob(::Type{NFormula}, value::Any) = NFormula(1)

function Hom(value::Any, A::NFormula, B::NFormula)::Formulas
  nin, nout = A.n, B.n
  @assert nout == 1
  vars = gensyms(nin)
  term = if isa(value, Symbol) && nin >= 1
    Formula(value, vars...)
  else
    value
  end
  Formulas([term], vars)
end

gensyms(n::Int) = [ gensym() for i in 1:n ]
gensyms(n::Int, tag) = [ gensym(tag) for i in 1:n ]
gensyms(A::NFormula, args...) = gensyms(A.n, args...)

""" Simultaneous substitution of variables in formula.
"""
function substitute(form::Formula, subst::Dict)
  Formula(head(form), [substitute(arg, subst) for arg in args(form)]...)
end

""" Simultaneous substitution of symbols in Julia expression.
"""
function substitute(expr::Expr, subst::Dict)
  Expr(expr.head, [substitute(arg, subst) for arg in expr.args]...)
end
substitute(sym::Symbol, subst::Dict) = get(subst, sym, sym)
substitute(x::Any, subst::Dict) = x

""" Create formula for sum of zero, one, or more terms.
"""
function sum_formulas(T::Type, terms::Vector)
  if length(terms) == 0
    zero(T)
  elseif length(terms) == 1
    terms[1]
  else
    Formula(:+, terms...)
  end
end
sum_formulas(terms::Vector) = sum_formulas(Int, terms)

""" Convert a formula, or formulas, to a wiring diagram.

The wiring diagram has an input for each symbol in `vars` and an output for
each given formula. All terminal symbols not appearing in `vars` are treated as
symbolic constants.

The algorithm creates wiring diagrams in normal form for a cartesian category,
meaning that subformulas are maximally shared (cf. `normalize_copy!` for wiring
diagrams and the congruence closure algorithm).
"""
function to_wiring_diagram(formula::Formula, vars::Vector{Symbol}; kw...)
  to_wiring_diagram([formula], vars; kw...)
end

function to_wiring_diagram(formulas::Vector{Formula}, vars::Vector{Symbol};
                           to_box::Function=to_untyped_box,
                           to_port::Function=to_untyped_port)
  diagram = WiringDiagram(map(to_port, vars), map(to_port, formulas))

  term_map = Dict{Union{Formula,Symbol},Tuple{Int,Int}}((
    var => (input_id(diagram), i) for (i, var) in enumerate(vars)
  ))
  function visit(term::Formula)
    get!(term_map, term) do
      preds = map(visit, args(term))
      v = add_box!(diagram, to_box(term))
      add_wires!(diagram, (pred => (v, i) for (i, pred) in enumerate(preds)))
      (v,1)
    end
  end
  function visit(sym::Symbol)
    get!(term_map, sym) do
      v = add_box!(diagram, to_box(sym))
      (v,1)
    end
  end
  function visit(num::Number)
    v = add_box!(diagram, to_box(num))
    (v,1)
  end

  outs = map(visit, formulas)
  add_wires!(diagram,
    (out => (output_id(diagram), i) for (i, out) in enumerate(outs)))
  diagram
end

to_untyped_box(form::Formula) =
  Box(head(form), repeat([nothing], length(args(form))), [nothing])
to_untyped_box(x) = Box(x, Nothing[], [nothing])
to_untyped_port(x) = nothing

# Evaluation
############

""" Evaluate a formula, optionally with vectorization.
"""
function evaluate(form::Formula, env::Dict; vector::Bool=true)
  evaluate_formula(form, env; vector=vector)
end

function evaluate_formula(form::Formula, env::Dict; vector::Bool=true)
  f = evaluate_formula(head(form), env)
  f_args = (evaluate_formula(arg, env) for arg in args(form))
  vector ? f.(f_args...) : f(f_args...)
end

function evaluate_formula(name::Symbol, env::Dict; kw...)
  # Look up name in formula environment.
  get(env, name) do
    # Failing that, look up name as field of Julia module.
    # XXX: Should the module always be Main?
    getfield(Main, name)
  end
end

evaluate_formula(literal::Number, env; kw...) = literal

# Pretty-print
##############

function show(io::IO, form::Formula)
  print(io, "Formula(")
  join(io, [ sprint(show, x) for x in [head(form); args(form)] ], ", ")
  print(io, ")")
end

function show(io::IO, ::MIME"text/latex", form::Formula)
  print(io, "\$")
  show_latex(io, form)
  print(io, "\$")
end

""" Show the expression in infix notation using LaTeX math.

Does *not* include `\$` or `\\[begin|end]{equation}` delimiters.
"""
show_latex(form::Formula) = show_latex(stdout, form)
show_latex(io::IO, form::Formula) = show_latex_formula(io, form)

""" Show the formula as an S-expression.

Cf. the standard library function `Meta.show_sexpr`.
"""
show_sexpr(form::Formula) = show_sexpr(stdout, form)

function show_sexpr(io::IO, form::Formula)
  print(io, "(")
  join(io, [sprint(show, head(form));
            [sprint(show_sexpr, arg) for arg in args(form)]], " ")
  print(io, ")")
end

# Show terminal nodes as LaTeX.

function show_latex_formula(io::IO, bool::Bool; kw...)
  print(io, bool ? "\\top" : "\\bot")
end

function show_latex_formula(io::IO, num::Number; kw...)
  if num < 0
    # Treat negative literals as compound expressions.
    show_latex_formula(io, Formula(:-, abs(num)); kw...)
  else
    print(io, isinf(num) ? "\\infty" : num)
  end
end

function show_latex_formula(io::IO, sym::Symbol; kw...)
  if sym in latex_symbol_aliases
    print(io, "\\", sym)
  else
    print(io, get(latex_symbol_table, sym) do
      length(string(sym)) == 1 ? sym : "\\mathrm{$sym}"
    end)
  end
end

# Show non-terminal nodes as LaTeX.

function show_latex_formula(io::IO, form::Formula; kw...)
  # Special case: Dispatch on special LaTeX commands (e.g., \frac).
  headsym = head(form)
  if haskey(latex_command_table, headsym)
    return show_latex_command(io, form, headsym; kw...)
  end

  # Special case: Dispatch on operators (prefix, infix, and postfix).
  nargs = length(args(form))
  if nargs == 1
    if haskey(latex_prefix_table, headsym)
      return show_latex_prefix(io, form, headsym; kw...)
    end
    if haskey(latex_postfix_table, headsym)
      return show_latex_postfix(io, form, headsym; kw...)
    end
  elseif nargs > 1
    if haskey(latex_infix_table, headsym)
      return show_latex_infix(io, form, headsym; kw...)
    end
  end

  # General case: Dispatch on head symbol as value type.
  show_latex_formula(io, form, Val{headsym}; kw...)
end

# Default LaTeX pretty-printer.

function show_latex_formula(io::IO, form::Formula, ::Type; kw...)
  show_latex_formula(io, head(form))
  print(io, "\\left(")
  join(io, [sprint(show_latex_formula, arg) for arg in args(form)], ",")
  print(io, "\\right)")
end

# Special LaTeX pretty-printers.

function show_latex_formula(io::IO, form::Formula, ::Type{Val{:^}};
                            parent_precedence::Int=typemax(Int), kw...)
  @assert length(args(form)) == 2
  precedence = 2
  latex_paren_when(io, precedence >= parent_precedence) do
    show_latex_formula(io, first(form); parent_precedence=precedence, kw...)
    print(io, "^{")
    show_latex_formula(io, last(form))
    print(io, "}")
  end
end
function show_latex_formula(io::IO, form::Formula, ::Type{Val{:.^}}; kw...)
  show_latex_formula(io, form, Val{:^}; kw...)
end

# LaTeX pretty-printers for unary and binary operators.

function show_latex_infix(io::IO, form::Formula, sym::Symbol;
                          parent_precedence::Int=typemax(Int), kw...)
  op, precedence = latex_infix_table[sym]
  sep = if isempty(op)
      any(isa(x,Number) for x in args(form)[2:end]) ? " \\cdot " : " "
    else
      " $op "
  end
  show_child = (io, form) ->
    show_latex_formula(io, form; parent_precedence=precedence, kw...)
  latex_paren_when(io, precedence >= parent_precedence) do
    join(io, [sprint(show_child, arg) for arg in args(form)], sep)
  end
end

function show_latex_prefix(io::IO, form::Formula, sym::Symbol;
                           parent_precedence::Int=typemax(Int), kw...)
  @assert length(args(form)) == 1
  op, precedence = latex_prefix_table[sym]
  latex_paren_when(io, precedence >= parent_precedence) do
    print(io, op, " ")
    show_latex_formula(io, first(form); parent_precedence=precedence, kw...)
  end
end

function show_latex_postfix(io::IO, form::Formula, sym::Symbol;
                            parent_precedence::Int=typemax(Int), kw...)
  @assert length(args(form)) == 1
  op, precedence = latex_postfix_table[sym]
  latex_paren_when(io, precedence >= parent_precedence) do
    show_latex_formula(io, first(form); parent_precedence=precedence, kw...)
    print(io, " ", op)
  end
end

function show_latex_command(io::IO, form::Formula, sym::Symbol; kw...)
  cmd = latex_command_table[sym]
  print(io, cmd, "{")
  join(io, [sprint(show_latex_formula, arg) for arg in args(form)], "}{")
  print(io, "}")
end

function latex_paren_when(f::Function, io::IO, paren::Bool)
  if (paren) print(io, "\\left(") end
  f()
  if (paren) print(io, "\\right)") end
end

# Operator precedence follows the usual conventions of mathematics.
#
# https://en.wikipedia.org/wiki/Order_of_operations
# https://rosettacode.org/wiki/Operator_precedence#Mathematica

const latex_infix_table = Dict{Symbol,Tuple{String,Int}}(
  #:^ => ("^", 2),
  #:.^ => ("^", 2),
  :* => ("", 4),
  :.* => ("", 4),
  :+ => ("+", 5),
  :- => ("-", 5),
  :(==) => ("=", 6),
  :!= => ("\\neq", 6),
  :< => ("<", 6),
  :> => (">", 6),
  :<= => ("\\leq", 6),
  :>= => ("\\geq", 6),
  :in => ("\\in", 6),
  :isa => (":", 6),
  :& => ("\\wedge", 8),
  :| => ("\\vee", 8),
  :-> => ("\\mapsto", 9),
)
const latex_prefix_table = Dict{Symbol,Tuple{String,Int}}(
  :- => ("-", 3),
  :! => ("\\neg", 3),
)
const latex_postfix_table = Dict{Symbol,Tuple{String,Int}}(
  :factorial => ("!", 2),
)

const latex_command_table = Dict{Symbol,String}(
  :/ => "\\frac",
  :./ => "\\frac",
  :binomial => "\\binom",
  :sqrt => "\\sqrt",
)
const latex_symbol_table = Dict{Symbol,String}(
  :Inf => "\\infty",
  :Int => "\\mathbb{Z}", :Integer => "\\mathbb{Z}",
  :Signed => "\\mmathbb{Z}", :Unsigned => "\\mathbb{N}",
  :Real => "\\mathbb{R}", :Complex => "\\mathbb{C}",
)
const latex_symbol_aliases = Set{Symbol}([
  :alpha, :beta, :gamma, :Gamma, :delta, :Delta, :epsilon, :varepsilon,
  :zeta, :eta, :theta, :vartheta, :Theta, :iota, :kappa, :lambda, :Lambda,
  :mu, :nu, :xi, :Xi, :pi, :Pi, :rho, :varrho, :sigma, :Sigma, :tau,
  :upsilon, :Upsilon, :phi, :varphi, :Phi, :chi, :psi, :Psi, :omega, :Omega,
  :exp, :log, :sin, :cos, :tan,
])

end


text:
using Preprocessing

function createcggrid(dgnodes,tlocal)

# Get some dimensions of the mesh
~, nd, ne = size(dgnodes);
~, nve    = size(tlocal);

cgelcon,~,~,~,cgnodes = Preprocessing.mkcgent2dgent(dgnodes,1e-7);

cgnodes = cgnodes';
cgcells = createcgcells(cgelcon,tlocal,ne);

elemtype = 0;
if (nd==2) && (nve==4)
    elemtype=1;
end
if (nd==3) && (nve==8)
    elemtype=1;
end

cell_t = getcelltype(nd,elemtype);

return cgnodes, cgelcon, cgcells, cell_t

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
using DataStructures

nested_dict(K::DataType=Any, V::DataType=Any) = DefaultDict{K, V}(nested_dict)

"""
    Path(s::String)

Create a path from a `String`. Paths may be concatenated with `*`
"""
struct Path
    x::String
end

Base.convert(::Type{String}, p::Path) = p.x
Base.String(p::Path) = convert(String, p)
Base.open(x::Path, m::String) = open(convert(String, x), m)

# Necessary to specify Function due to ambiguity errors
Base.open(f::Function, x::Path) = open(f, convert(String, x))
Base.:*(x::Path, y::Union{Path, String}) = Path(joinpath(String(x), String(y)))
Base.split(x::Path) = splitpath(x.x)

for f in [:(Base.open), :(Base.splitpath)]
    @eval $f(x::Path) = $f(String(x))
end

for f in [:(Base.basename), :(Base.isdir), :(Base.ispath), :(Base.isfile)]
    @eval $f(x::Path) = Path($f(String(x)))
end

"""
    @indir directory expression

Move to `directory`, execute expression, then move back.
"""
macro indir(directory, ex)
    quote
        local dir = pwd()
        cd($(directory))
        local val = $(esc(ex))
        cd(dir)
        val
    end
end

max_stars_repo_name: kescobo/BadBio

max_issues_repo_name: kescobo/BadBio

max_forks_repo_name: kescobo/BadBio


text:
export simdir, simid, @run, @runsync, @rerun, @rerunsync, in_simulation_mode, @SimulationEnvironment

const ENV_SIM_FOLDER = "SIMULATION_FOLDER"
const ENV_SIM_ID = "SIMULATION_ID"

abstract type AbstractSimulationEnvironment end

struct DefaultSimulation <: AbstractSimulationEnvironment end

macro SimulationEnvironment(name::Symbol)
    :(struct $name <: DrWatsonSim.AbstractSimulationEnvironment end)
end

from_folder_name(n::String) = parse(Int, n)
to_folder_name(n) = string(n)
in_simulation_mode() = ENV_SIM_ID in keys(ENV)

function simdir(args...)
    if ENV_SIM_FOLDER in keys(ENV)
        return joinpath(ENV[ENV_SIM_FOLDER],args...)
    end
    error("Not in simulation environment")
end

function simid()
    if ENV_SIM_ID in keys(ENV)
        return parse(Int,ENV[ENV_SIM_ID])
    end
    error("Not in simulation environment")
end

function get_next_simulation_id(folder)
    id = 1
    for i in 1:100_000_000
        try
            mkdir(joinpath(folder,to_folder_name(id)))
            return id
        catch e
            if !isa(e, Base.IOError)
                rethrow(e)
            end
            id += 1
        end
    end
    error("Couldn't genereate new id in '$folder'")
end

run_simulation(f,p,args...) = run_simulation(f, [p], args...)

function run_in_simulation_mode(f)
    m = Metadata(simdir())
    f(m["parameters"])
end

function add_simulation_metadata!(::AbstractSimulationEnvironment,
                                  m::Metadata;
                                 simulation_id,
                                 simulation_submit_time,
                                 simulation_submit_group,
                                 scriptfile,
                                 command,
                                 env,
                                 source
                                )
    m["simulation_submit_time"] = simulation_submit_time
    m["simulation_submit_group"] = simulation_submit_group
    m["simulation_id"] = simulation_id
    m["mtime_scriptfile"] = mtime(scriptfile)
    m["julia_command"] = command
    m["ENV"] = env
    tag!(m, source=source)
end

submit_command(::AbstractSimulationEnvironment,id,env) = `$(Base.julia_cmd()) $(PROGRAM_FILE)`

function run_simulation(t::AbstractSimulationEnvironment,f,param,directory,source; wait_for_finish=false)
    if in_simulation_mode()
        run_in_simulation_mode(f)
    else
        simulation_ids = map(param) do _
            get_next_simulation_id(directory)
        end

        simulation_submit_time = DrWatsonSim.now()
        simulation_submit_group = [standardize_path(joinpath(directory,to_folder_name(i))) for i in simulation_ids]

        tasks = map(zip(param,simulation_ids)) do (p,id)
            folder = joinpath(directory,to_folder_name(id))
            m = Metadata(folder)
            env = copy(ENV)

            env[ENV_SIM_FOLDER] = folder
            env[ENV_SIM_ID] = string(id)

            command = submit_command(t, id, env)

            # Add the metadata for the simulation mode.
            # This is the optional data.
            add_simulation_metadata!(t, m, simulation_id=id, simulation_submit_time=simulation_submit_time,
                                     simulation_submit_group=simulation_submit_group,
                                     scriptfile=PROGRAM_FILE, command=command, env=env, source=source)

            # Add the parameters as metadata.
            # This is the only metadata that is required
            m["parameters"] = p

            return @async begin
                run(setenv(`$command`, env), wait=wait_for_finish)
            end
        end
        print("Starting $(length(simulation_ids)) job(s):")
        for id in simulation_ids
            project_rel_path = standardize_path(joinpath(directory,to_folder_name(id)))
            print("\n  $project_rel_path")
        end
        println()
        wait.(tasks)
    end
end

function rerun_simulation(t::AbstractSimulationEnvironment,f,folder,source; wait_for_finish=false)
    if in_simulation_mode()
        run_in_simulation_mode(f)
    else
        m_original = Metadata(folder)
        simulation_submit_time = DrWatsonSim.now()
        simulation_submit_group = "simulation_submit_group" in keys(m_original) ? m_original["simulation_submit_group"] : []
        p = m_original["parameters"]
        id = from_folder_name(basename(folder))
        m = Metadata!(folder)

        env = copy(ENV)

        env[ENV_SIM_FOLDER] = folder
        env[ENV_SIM_ID] = string(id)

        command = submit_command(t, id, env)

        add_simulation_metadata!(t, m, simulation_id=id, simulation_submit_time=simulation_submit_time,
                                 simulation_submit_group=simulation_submit_group,
                                 scriptfile=PROGRAM_FILE, command=command, env=env, source=source)

        m["parameters"] = p
        t = @async run(setenv(`$command`, env), wait=wait_for_finish)
        wait(t)
    end
end

macro run(f, p, directory)
    source=QuoteNode(__source__)
    :(run_simulation(DrWatsonSim.DefaultSimulation(),$(esc(f)), $(esc(p)), $(esc(directory)), $source))
end

macro runsync(f, p, directory)
    source=QuoteNode(__source__)
    :(run_simulation(DrWatsonSim.DefaultSimulation(),$(esc(f)), $(esc(p)), $(esc(directory)), $source, wait_for_finish=true))
end

macro rerun(f, path)
    source=QuoteNode(__source__)
    :(rerun_simulation(DrWatsonSim.DefaultSimulation(),$(esc(f)), $(esc(path)), $source))
end

macro rerunsync(f, path)
    source=QuoteNode(__source__)
    :(rerun_simulation(DrWatsonSim.DefaultSimulation(),$(esc(f)), $(esc(path)), $source, wait_for_finish=true))
end

macro run(t, f, p, directory)
    source=QuoteNode(__source__)
    :(run_simulation($(esc(t)),$(esc(f)), $(esc(p)), $(esc(directory)), $source))
end

macro runsync(t, f, p, directory)
    source=QuoteNode(__source__)
    :(run_simulation($(esc(t)),$(esc(f)), $(esc(p)), $(esc(directory)), $source, wait_for_finish=true))
end

macro rerun(t, f, path)
    source=QuoteNode(__source__)
    :(rerun_simulation($(esc(t)),$(esc(f)), $(esc(path)), $source))
end

macro rerunsync(t, f, path)
    source=QuoteNode(__source__)
    :(rerun_simulation($(esc(t)),$(esc(f)), $(esc(path)), $source, wait_for_finish=true))
end


text:
import FractionalTimeDG.P, FractionalTimeDG.dP

α = 3/4
r = 4
store = FractionalTimeDG.Store(α, r, 2r)

function brute_force_Hn(r::Integer, n::Integer, ℓbar::Integer, α::T,
                        t::OffsetVector{T},
                        rtol=1e-8) where T <: AbstractFloat

    function inner(τ::T, j::Integer)
        I, err = quadgk(-one(T), one(T), rtol=rtol) do σ
            return ( 1 + Δ(n, ℓbar, τ, σ, t) )^(α-1) * P(j-1, σ)
        end
        return I
    end

    Hnℓbar = Array{T}(undef, r, r)
    ℓ = n - ℓbar
    kℓ = t[ℓ] - t[ℓ-1]
    D = (t[n]+t[n-1])/2 - (t[ℓ]+t[ℓ-1])/2
    c = D^(α-1) * kℓ / ( 2Γ(α) )
    for j = 1:r
        Aj = inner( one(T), j)
        Bj = inner(-one(T), j)
        pow = one(T)
        for i = 1:r
            Cij, err = quadgk(-one(T), one(T), rtol=rtol) do τ
                dP(i-1,τ) * inner(τ, j)
            end
            pow = -pow
            Hnℓbar[i,j] = c * ( Aj + pow * Bj - Cij )
        end
    end
    return Hnℓbar
end

N = 4
n = 3

t = OffsetVector([0.0, 0.35, 0.8, 1.3, 2.0], 0:N)

Hn = coef_Hn!(r, n, n-1, t, r+2, store)

bf_H1 = brute_force_Hn(r, n, 1, α, t)
err1 = bf_H1 - Hn[1]
@test all( abs.(err1) .< 1e-8 )

bf_H2 = brute_force_Hn(r, n, 2, α, t)
err2 = bf_H1 - Hn[2]
@test all( abs.(err1) .< 1e-8 )


text:
module IPythonDisplay

using IJulia
import IJulia: send_ipython, publish, msg_pub, execute_msg, display_dict, displayqueue, undisplay

import Base: display, redisplay
export display, redisplay, InlineDisplay, undisplay

immutable InlineDisplay <: Display end

# supported MIME types for inline display in IPython, in descending order
# of preference (descending "richness")
const ipy_mime = [ "text/html", "text/latex", "image/svg+xml", "image/png", "image/jpeg", "text/plain" ]

for mime in ipy_mime
    @eval begin
        function display(d::InlineDisplay, ::MIME{symbol($mime)}, x)
            send_ipython(publish,
                         msg_pub(execute_msg, "display_data",
                                 ["source" => "julia", # optional
                                  "metadata" => Dict(), # optional
                                  "data" => [$mime => stringmime(MIME($mime), x)] ]))
        end
    end
end

# deal with annoying application/x-latex == text/latex synonyms
display(d::InlineDisplay, m::MIME"application/x-latex", x) = display(d, MIME("text/latex"), stringmime(m, x))

# override display to send IPython a dictionary of all supported
# output types, so that IPython can choose what to display.
function display(d::InlineDisplay, x)
    undisplay(x) # dequeue previous redisplay(x)
    send_ipython(publish,
                 msg_pub(execute_msg, "display_data",
                         ["source" => "julia", # optional
                          "metadata" => Dict(), # optional
                          "data" => display_dict(x) ]))
end

# we overload redisplay(d, x) to add x to a queue of objects to display,
# with the actual display occuring when display() is called or when
# an input cell has finished executing.

function redisplay(d::InlineDisplay, x)
    if !contains(displayqueue, x)
        push!(displayqueue, x)
    end
end

function display()
    q = copy(displayqueue)
    empty!(displayqueue) # so that undisplay in display(x) is no-op
    for x in q
        display(x)
    end
end

end # module


text:
# # Burgers (Liquid)
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/model_burgers.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/model_burgers.ipynb)

using RHEOS
## include a helper function for plotting
include("assets/plothelper.jl");
#-

BurgersLiquid

# #### Constitutive Equation
# ```math
# \sigma(t) + p_1 \frac{d \sigma(t)}{dt} + p_2  \frac{d^2 \sigma(t)}{dt^2} = q_1  \frac{d \epsilon(t)}{dt} + q_2  \frac{d^2 \epsilon(t)}{dt^2}
# ```

# ```math
# \text{where}\; \ p_1 = \frac{\eta_1}{k_1}+\frac{\eta_1}{k_2}+\frac{\eta_1}{k_2}\text{,}\; \ p_2 = \frac{\eta_1 \eta_2}{k_1 k_2}\text{,}\; \ q_1 =  \eta_1\; \ \text{and}\; \ q_2 = \frac{\eta_1 \eta_2}{k_2}
# ```

# #### Relaxation Modulus
# ```math
# G(t) = \left[ \left(q_1 - q_2 r_1\right) e^{-r_1 t}  - \left(q_1 - q_2 r_2\right) e^{-r_2 t} \right]/A
# ```

# ```math
# \text{where}\; \ r_1 = (p_1-A)/2p_2 \text{,}\; \ r_2 = (p_1 +A) /2p_2 \text{, and}\; \ A = \sqrt{p_1^2 - 4p_2}
# ```

# #### Creep Modulus
# ```math
# J(t) = \frac{1}{k_1} + \frac{t}{\eta_1} + \frac{1}{k_2} \left(1-e^{-k_2 t/\eta_2}\right)
# ```

# #### Storage Modulus
# ```math
# G^{\prime}(\omega) = \frac{p_1 q_1 \omega^2 - q_2 \omega^2 (1-p_2 \omega^2)}{p_1^2 \omega^2 + (1-p_2 \omega^2)^2}
# ```

# #### Loss Modulus
# ```math
# G^{\prime \prime}(\omega) = \frac{p_1 q_2 \omega^3 + q_1 \omega (1-p_2 \omega^2)}{p_1^2 \omega^2 + (1-p_2 \omega^2)^2}
# ```

# ## Qualitative Behaviours of the Moduli

models = Vector{RheoModel}()

## plot moduli for varying k₁
for k in [0.5, 1.0, 2.0]

    push!(models, RheoModel(BurgersLiquid, (η₁ = 1, k₁ = k, η₂ = 1, k₂ = 1)))

end

plotmodel(models)

# **Reference**: Findley, William N., and Francis A. Davis. Creep and relaxation of nonlinear viscoelastic materials. Courier Corporation, 2013.


text:
using DelimitedFiles, Econometrics
function main()
# Check the Nerlove data for collinearity
data = readdlm("nerlove.data")
data = data[:,3:6] # these are the regressors
data = log.(data)
n = size(data,1)

# Indices of remaining regressors, after selecting each in turn
j = [2 3 4; 1 3 4; 1 2 4; 1 2 3]
title = ["output", "labor","capital","fuel"]
for i = 1:4
	println()
	println(title[i])
	y = data[:,i] # select the artificial dep. var.
	k = j[i,:] # select the other regressors
	names = ["constant"; title[k]]
	x = [ones(n,1) data[:,k]]
	ols(y, x, names=names)
end
return
end
main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
import JSON

## As simple as...
## (run netexp)
## formatbioxpoutput(netexp_simple_results_dir)
## OR
## formatbioxpoutput(netexp_simple_results_file)
"""
Convert JGI files from old sytle (used in original Seeding life submission) to
new style (used for BioXP rebuild)
"""
function old_to_new_rids(olddir::String,newdir::String)
    for f in readdir(olddir)

        oldjson = JSON.parsefile(joinpath(olddir,f))
        ## Make sure keys in reactions and products are the same
        @assert Set(keys(oldjson["substrates"])) == Set(keys(oldjson["products"]))

        if !ispath(newdir)
            mkpath(newdir)
        end

        newpath = joinpath(newdir,f)
        open(newpath,"w") do f
            JSON.print(f,keys(oldjson["substrates"]),2)
        end

    end
end

## Single file
function formatbioxpoutput(path::String,
                           write_path::Union{String,Nothing}=nothing)

    if isdir(path)
        _formatbioxpdir(path,write_path)
    elseif isfile(path)
        _formatbioxpfile(path,write_path)
    else
        error("neither a file nor a directory")
    end

end

function _formatbioxpfile(fname::String,
                         write_path::Union{String,Nothing}=nothing)

    ## write_path is guaranteed to be a file here, named "formatted/fname.ext"
    if isnothing(write_path)
        # path = joinpath(splitpath(abspath(path))[1:end-2]...)
        newdir = joinpath(dirname(abspath(fname)),"formatted")
        if !ispath(newdir)
            mkpath(newdir)
        end
        write_path = joinpath(newdir,basename(fname))
    ## in case write_path is a file, get the directory
    elseif !ispath(dirname(abspath(write_path)))
        mkpath(dirname(abspath(write_path)))
    end

    formattedexpansion = formatexpansion(readexpansion(fname))

    open(write_path,"w") do f
        JSON.print(f, formattedexpansion, 2) #indent=2
    end
end

function _formatbioxpdir(dir::String,
                        write_path::Union{String,Nothing}=nothing)

    ## write_path is guaranteed to be a dir here, named "formatted"
    if isnothing(write_path)
        newdir = joinpath(abspath(dir),"formatted")
        if !ispath(newdir)
            mkpath(newdir)
        end
    ## This can only possibly create a directory
    elseif !ispath(write_path)
        mkpath(write_path)
    end

    # "$i.json"
    for fname in filter(x->endswith(x,".json"),readdir(abspath(dir))) ## List all json files in dir

        formattedexpansion = formatexpansion(readexpansion(joinpath(abspath(dir),fname)))

        if isnothing(write_path)
            real_write_path = joinpath(newdir,basename(fname))
        ## Only supports choosing the dirname to write when writing whole dir
        else
            real_write_path = joinpath(write_path,basename(fname))
        end

        open(real_write_path,"w") do f
            JSON.print(f, formattedexpansion, 2) #indent=2
        end
    end

end

function formatexpansion(nexp::Expansion)
    formatexpansion(nexp.x,
                    nexp.t,
                    nexp.compounds,
                    nexp.reactions,
                    nexp.X,
                    nexp.Y)
end

function formatexpansion(x::Vector{Int},t::Vector{Int},compounds::Array{String,1},reactions::Array{String,1},X::Vector{Vector{Int}},Y::Vector{Vector{Int}})
    data = Dict("stats"=>Dict(),"generations"=>Dict())

    ## Store data about the scope, about the equilibrium
    data["stats"]["scope_compounds"] = compounds
    data["stats"]["scope_reactions"] = reactions
    data["stats"]["scope_seeds"] = [compounds[i] for i in 1:length(x) if Bool(x[i])]
    data["stats"]["scope_targets"] = [compounds[i] for i in 1:length(t) if Bool(t[i])]

    # [Int(i in seed_list) for i in compounds]

    ## Store data by generation
    n_generations = length(X)-1 ## Last two generations are repeated for compounds
    for gen in 1:n_generations

        data["generations"][gen] = Dict()
        # data["generations"][gen]["reactions_cumulative"] = []
        data["generations"][gen]["compounds_cumulative"] = []
        data["generations"][gen]["targets_cumulative"] = []
        data["generations"][gen]["reactions_new"] = []
        data["generations"][gen]["compounds_new"] = []
        data["generations"][gen]["targets_new"] = []

        ## Add cumulative reaction list for each generation
        data["generations"][gen]["reactions_cumulative"] = [reactions[i] for i in 1:length(reactions) if Bool(Y[gen][i])]
        ## Alternative way to write
        # for (i,r) in enumerate(reactions)
        #     if Y[gen][i]==1
        #         push!(data["generations"][gen]["reactions_cumulative"],reactions[i])
        #     end
        # end

        ## Add cumulative compound and target list for each generation
        # data["generations"][gen]["compounds_cumulative"] = [compounds[i] for i in 1:length(compounds) if Bool(X[gen][i])]
        for (i,c) in enumerate(compounds)
            if X[gen][i]==1
                push!(data["generations"][gen]["compounds_cumulative"],compounds[i])
                if t[i]==1
                    push!(data["generations"][gen]["targets_cumulative"],compounds[i])
                end
            end
        end

        ## Store new reactions and compounds
        if gen==1
            data["generations"][gen]["reactions_new"] = collect(Set(data["generations"][gen]["reactions_cumulative"]))
            data["generations"][gen]["compounds_new"] = collect(Set(data["generations"][gen]["compounds_cumulative"]))
            data["generations"][gen]["targets_new"] = collect(Set(data["generations"][gen]["targets_cumulative"]))
        end

        if gen!=1
            data["generations"][gen]["reactions_new"] = setdiff(data["generations"][gen]["reactions_cumulative"], data["generations"][gen-1]["reactions_cumulative"])
            data["generations"][gen]["compounds_new"] = setdiff(data["generations"][gen]["compounds_cumulative"], data["generations"][gen-1]["compounds_cumulative"])
            data["generations"][gen]["targets_new"] = setdiff(data["generations"][gen]["targets_cumulative"], data["generations"][gen-1]["targets_cumulative"])
        end
    end

    data
end

max_stars_repo_name: hbsmith/BioXP

max_issues_repo_name: hbsmith/BioXP

max_forks_repo_name: hbsmith/BioXP


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

using Markdown
using InteractiveUtils

# ╔═╡ cf9eb2f0-1eb9-11ec-2557-37ce61a5e120
using RemoteS, GMT

# Interpolate from sensor to geographic coordinates at approx 1 km grid steps
G = grid_at_sensor("C:/v/AQUA_MODIS.20210805T131001.L2.SST.NRT.nc", "sst", inc=0.01);

# ╔═╡ 8bc57780-1eb9-11ec-0581-7b462f0ed4eb
md"# Interpolate a MODIS L2 SST file

The Level 2 MODIS files are the ones that contain the information at its maximum spatial resolution. But in those products the data is not equally spaced so we must interpolate it first before use. The `grid_at_sensor` function provides a handy interface to do it easily.

Imagine that you have downloaded the file `AQUA_MODIS.20210805T131001.L2.SST.NRT.nc` from the [OceanColor](https://oceancolor.gsfc.nasa.gov/) site (eventually following the example in [Plot AQUA satellite tracks](https://www.generic-mapping-tools.org/RemoteS.jl/dev/gallery/Aqua_orbits/remotes_sat_tracks/))
"

# ╔═╡ 13755650-1eba-11ec-38bc-ed6e4c29d316
# Display it
imshow(G, proj=:guess, coast=true, dpi=200)

# ╔═╡ 41c73f50-1eba-11ec-3c28-1fc77965fbda
md"Note that since we used a grid step 0f 0.01 for the interpolation and this value is very close to the MODIS maximum spatial resolution, the left and right regions have beam spacing 2 to 5 times this and show many little holes. Do not confuse these little holes with the larger ones that are caused by cloud coverage. So we will recalculate the grid at increments of ~2 km and over the region that has the higher data contents.
"

# ╔═╡ 19448a90-1ebc-11ec-1b1e-53bcc9df7538
# Recompue at a inc=0.02 and over a sub-region
G = grid_at_sensor("C:/v/AQUA_MODIS.20210805T131001.L2.SST.NRT.nc", "sst", region=(-13,10,33.8,44.5), inc=0.02);

# Make a nicer image with illumination.
imshow(G, proj=:guess, coast=true, shade=true, title="Sea Surace Temperature", colorbar=true)

# ╔═╡ Cell order:
# ╟─8bc57780-1eb9-11ec-0581-7b462f0ed4eb
# ╠═cf9eb2f0-1eb9-11ec-2557-37ce61a5e120
# ╠═13755650-1eba-11ec-38bc-ed6e4c29d316
# ╠═41c73f50-1eba-11ec-3c28-1fc77965fbda
# ╠═19448a90-1ebc-11ec-1b1e-53bcc9df7538

max_stars_repo_name: GenericMappingTools/RemoteS

max_issues_repo_name: GenericMappingTools/RemoteS

max_forks_repo_name: GenericMappingTools/RemoteS


text:
include("AWSCorePrototypeServices.jl")
s3 = AWSCorePrototypeServices.s3

function low_level_main()
    buckets = s3("GET", "/")
    println(buckets)
end

max_stars_repo_name: nicoleepp/JuliaCloud-AWS-State

max_issues_repo_name: nicoleepp/JuliaCloud-AWS-State

max_forks_repo_name: nicoleepp/JuliaCloud-AWS-State


text:
module Synchronous

using ..Instantiation

import ModiaMath
import ModiaMath.ModiaToModiaMath  # ModiaSimulationModel, EventRestart, Restart, positive!, ...
import ModiaMath.ModiaToModiaMath.ModiaSimulationModel
import ModiaMath.Restart

@static if VERSION < v"0.7.0-DEV.2005"
    Nothing = Void
end

#= Synchronous Modelica primitives

Clock()
Clock(intervalCounter::Int, resolution::Int=1)  # Add version with rational number
Clock(interval::Float64)
Clock(condition::Bool, startInterval::Float64=0.0)
Clock(c::Clock, solverMethod)

previous(x)

sample(u, c::Clock)
sample(u)

hold(u) = u

subSample(u, factor::Int)
superSample(u, factor::Int)
shiftSample(u, shiftCounter::Int, resolution::Int=1)  # Add version with rational number
backSample(u, backCounter::Int, resolution::Int=1)  # Add version with rational number

noClock()

interval(u)

=#

# The argument list of function F is extended with the argument simulationModel.

const nClock = 100
const nSample = 100
const nPrevious = 100
const nCrossings = 100

mutable struct Store
    clock::Vector{ModiaMath.Time}
    sample::Vector{Any}
    previous::Vector{Any}
    nextPrevious::Vector{Any}
    crossings::Vector{Any}

    @static if VERSION < v"0.7.0-DEV.2005"
      Store() = new(zeros(nClock), Array{Any}(nSample), zeros(nPrevious), zeros(nPrevious), zeros(nCrossings))
    else
      Store() = new(zeros(nClock), Array{Any}(undef, nSample), zeros(nPrevious), zeros(nPrevious), zeros(nCrossings))
    end
end

isLog(m) = true

# -------------------------------------------------------------------

# Temporary:
allInstances(s) = return s

function Clock(interval::ModiaMath.Time, m::ModiaSimulationModel, nr::Int)
    if ModiaToModiaMath.isPreInitial(m)
        if ModiaToModiaMath.isLogEvents(m)
            println("        in Clock, nr = $nr (isPreInitial)")
        end
        return false
    end

    if ModiaToModiaMath.isInitial(m)
        if ModiaToModiaMath.isLogEvents(m)
            println("        in Clock, nr = $nr (isInitial)")
        end

        if typeof(m.store) == Nothing
            m.store = Store()
        end

        m.store.clock[nr] = ModiaMath.getTime(m)
        ModiaToModiaMath.setNextEvent!(m, m.store.clock[nr])

    elseif ModiaToModiaMath.isEvent(m) && ModiaToModiaMath.isAfterSimulationStart(m)
        tick = abs(ModiaMath.getTime(m) - m.store.clock[nr]) < 1E-10
        if tick
            if ModiaToModiaMath.isLogEvents(m)
                println("        in Clock, nr = $nr (isEvent; clock is active)")
            end

            m.store.clock[nr] = ModiaMath.getTime(m) + interval
            ModiaToModiaMath.setNextEvent!(m, m.store.clock[nr])
            return true
        end
        ModiaToModiaMath.setNextEvent!(m, m.store.clock[nr])
    end

    return false
end

function positive(crossing::Float64, m::ModiaSimulationModel, nr::Int; restart::ModiaToModiaMath.EventRestart=Restart)
    if typeof(m.store) == Nothing
        m.store = Store()
    end
        println("positive:"); @show nr m; ModiaToModiaMath.positive!(crossing, m, nr; restart=restart)
end

positiveChange(crossing::Float64, m::ModiaSimulationModel, nr::Int; restart::ModiaToModiaMath.EventRestart=Restart) =
         ModiaToModiaMath.positiveChange!(crossing, m, nr; restart=restart)
positiveEdge(crossing::Float64, m::ModiaSimulationModel, nr::Int; restart::ModiaToModiaMath.EventRestart=Restart) =
         ModiaToModiaMath.positiveEdge!(crossing, m, nr; restart=restart)

function sample(v, clock::Bool, m::ModiaSimulationModel, nr::Int)
    if ModiaToModiaMath.isPreInitial(m)
        if ModiaToModiaMath.isLogEvents(m)
            println("        in sample, nr = $nr (isPreInitial)")
        end
        return v
    elseif ModiaToModiaMath.isInitial(m)
        if ModiaToModiaMath.isLogEvents(m)
            println("        in sample, nr = $nr (initialize sample store)")
        end

        if typeof(m.store) == Nothing
            m.store = Store()
        end

        m.store.sample[nr] = v
        return v
    end

    if clock
        if ModiaToModiaMath.isLogEvents(m)
            println("        in sample, nr = $nr (clock is active)")
        end
        m.store.sample[nr] = v;
        return v
    else
        return m.store.sample[nr]
    end
end

# initPrevious not used
function initPrevious(v, m::ModiaSimulationModel, nr::Int)
    if ModiaToModiaMath.isPreInitial(m)
        if ModiaToModiaMath.isLogEvents(m)
            println("        in initPrevious, nr = $nr (isPreInitial)")
        end
        return v

    elseif ModiaToModiaMath.isInitial(m)
        if ModiaToModiaMath.isLogEvents(m)
            println("        in initPrevious, nr = $nr (initialize previous store)")
        end

        if typeof(m.store) == Nothing
            m.store = Store()
        end

        m.store.previous[nr] = v
        return v
    else
        return m.store.previous[nr]
    end
end

function previous(v, clock::Bool, m::ModiaSimulationModel, nr::Int)
    if ModiaToModiaMath.isPreInitial(m)
        if ModiaToModiaMath.isLogEvents(m)
            println("        in previous, nr = $nr (isPreInitial)")
        end
        return v

    elseif ModiaToModiaMath.isInitial(m)
        if ModiaToModiaMath.isLogEvents(m)
            println("        in previous, nr = $nr (initialize previous store)")
        end

        if typeof(m.store) == Nothing
            m.store = Store()
        end

        m.store.previous[nr]     = v
        m.store.nextPrevious[nr] = v
        return v
    end

    if clock
        m.store.previous[nr] = m.store.nextPrevious[nr];
        if ModiaToModiaMath.isLogEvents(m)
            println("        in previous, nr = $nr (clock is active)")
        end
    end
    return m.store.previous[nr]
end

#=
function previous(v, clock::Bool)
  return v
end
=#

function updatePrevious(v, m::ModiaSimulationModel, nr::Int)
    if ModiaToModiaMath.isInitial(m)
        if typeof(m.store) == Nothing
            m.store = Store()
        end
    end

    if ModiaToModiaMath.isEvent(m) && ModiaToModiaMath.isAfterSimulationStart(m)
    # if ModiaToModiaMath.isLogEvents(m)
    #   println("        in updatePrevious, nr = $nr, v = $v, time = ", ModiaMath.getTime(m))
    # end
        m.store.nextPrevious[nr] = v
    end
end

hold(v) = v

end


text:
module Chains

using ArgCheck: @argcheck
using Parameters: @unpack
using DocStringExtensions: FUNCTIONNAME, SIGNATURES, TYPEDEF

####
#### Index layouts
####

"""
$(TYPEDEF)

Allows extraction of a scalar or an array with the given dimension from a vector, starting
at `offset + 1`.
"""
struct IndexLayout{D <: Tuple{Vararg{Int}}}
    offset::Int
    len::Int
    dims::D
    function IndexLayout(offset, len, dims::D) where D
        @argcheck offset ≥ 0
        @argcheck len > 0
        @argcheck prod(dims) == len
        new{D}(offset, len, dims)
    end
end

function Base.getindex(A::AbstractVector, layout::IndexLayout)
    @unpack offset, len, dims = layout
    if isempty(dims)
        A[offset + 1]
    else
        reshape(view(A, offset .+ (1:len)), dims)
    end
end

function _getindex_stage(A, row_index, col_index)
    ix, _ = to_indices(A, (row_index, Colon()))
    if ix isa Integer
        A[ix, :][col_index]
    else
        map(row -> row[col_index], eachrow(view(A, ix, :)))
    end
end

function Base.getindex(A::AbstractMatrix, row_index, layout::IndexLayout)
    _getindex_stage(A, row_index, layout)
end

# optimized case for scalar layouts
function Base.getindex(A::AbstractMatrix, row_index, layout::IndexLayout{Tuple{}})
    col_index = layout.offset + 1
    ix, _ = to_indices(A, (row_index, Colon()))
    if ix isa Integer
        A[ix, col_index]
    else
        view(A, ix, col_index)
    end
end

_dims_string(dims::Tuple{}) = "scalar"
_dims_string(dims::Tuple{Int}) = "Vector($(first(dims)))"
_dims_string(dims) = "Array$(dims)"

function Base.show(io::IO, layout::IndexLayout)
    @unpack offset, len, dims = layout
    printstyled(io, "$((offset+1):(offset+len)) "; color = :light_black)
    printstyled(io, _dims_string(dims); color = :cyan)
end

"""
Type alias for a NamedTuple of array dimensions. Not exported.
"""
const NamedDims = NamedTuple{N,T} where {N,T<:Tuple{Vararg{Tuple{Vararg{Int}}}}}

"""
$(SIGNATURES)

Calculate column layouts from a `NamedTuple`) of array sizes.
"""
function _column_layouts(named_dims::NamedDims)
    acc = 0
    layouts = map(named_dims) do dims
        len = prod(dims)
        layout = IndexLayout(acc, len, dims)
        acc += len
        layout
    end
    layouts
end

####
#### Column schema
####

"""
$(TYPEDEF)

A named, ordered collection of `IndexedLayout`s. `getindex` methods with a symbol use the
given layout, `:` extracts a `NamedTuple`.
"""
struct IndexSchema{T}
    layouts::T
    function IndexSchema(named_dims::NamedDims)
        layouts = _column_layouts(named_dims)
        new{typeof(layouts)}(layouts)
    end
end

function Base.show(io::IO, cs::IndexSchema)
    get(io, :typeinfo, true) && print(io, "Column schema with layouts")
    for (name, layout) in pairs(layouts(cs))
        printstyled(io, "\n    ", name; color = :blue)
        print(io, " ", layout)
    end
end

layouts(cs::IndexSchema) = getfield(cs, :layouts)

Base.getproperty(cs::IndexSchema, name::Symbol) = getfield(layouts(cs), name)

Base.propertynames(cs::IndexSchema) = propertynames(layouts(scs))

function Base.length(cs::IndexSchema)
    last_layout = layouts(cs)[end]
    last_layout.offset + last_layout.len
end

"""
$(SIGNATURES)

Return a vector of `(name, indexes...)` for each index in `IndexSchema`, where
`indexes == ()` for scalars.

Useful for going back from flat indexes to layouts.
"""
function labels(cs::IndexSchema)
    ret = Vector{Any}(undef, length(cs))
    for (name, layout) in pairs(layouts(cs))
        @unpack offset, dims = layout
        if isempty(dims)
            ret[offset + 1] = (name, )
        else
            for (i, ci) in enumerate(CartesianIndices(map(i -> 1:i, dims)))
                ret[offset + i] = (name, Tuple(ci)...)
            end
        end
    end
    ret
end

function Base.getindex(V::AbstractVector, cs::IndexSchema)
    @argcheck !Base.has_offset_axes(V)
    @argcheck length(cs) == length(V)
    map(layout -> V[layout], layouts(cs))
end

function Base.getindex(A::AbstractMatrix, row_index, cs::IndexSchema)
    _getindex_stage(A, row_index, cs)
end

####
#### Chain
####

struct Chain{S <: IndexSchema, M <: AbstractMatrix}
    schema::S
    sample_matrix::M
    thinning::Int
    warmup::Int
    is_ordered::Bool
    function Chain(schema::S, sample_matrix::M;
                   thinning::Int = 1, warmup::Int = 0,
                   is_ordered::Bool = false) where {S <: IndexSchema, M <: AbstractMatrix}
        @argcheck !Base.has_offset_axes(sample_matrix)
        @argcheck length(schema) == size(sample_matrix, 2)
        @argcheck thinning ≥ 1
        @argcheck warmup ≥ 0
        new{S,M}(schema, sample_matrix, thinning, warmup, is_ordered)
    end
end

"""
$(SIGNATURES)

A `NamedTuple` of `IndexLayout`s for views into the sample matrix.
"""
schema(chain::Chain) = chain.schema

labels(chain::Chain) = labels(schema(chain))

function Base.show(io::IO, chain::Chain)
    @unpack schema, sample_matrix, thinning, warmup, is_ordered = chain
    print(io, is_ordered ? "Ordered" : "Unordered", " MCMC chain of ")
    printstyled(io, "$(size(sample_matrix, 1)) rows"; color = :red)
    is_ordered && print(io, " of which $(warmup) are warmup, thinning $(thinning)")
    print(IOContext(io, :typeinfo => false), " with schema", schema)
end

function sample_matrix(chain::Chain, include_warmup::Val{false})
    @view chain.sample_matrix[(chain.warmup + 1):end, :]
end

sample_matrix(chain::Chain, include_warmup::Val{true}) = chain.sample_matrix

"""
$(FUNCTIONNAME)(chain, include_warmup = Val{false}())

The sample matrix, with or without the warmup.
"""
sample_matrix(chain::Chain) = sample_matrix(chain, Val{false}())

"""
$(SIGNATURES)

The thinning for the sample. When `is_ordered(chain) == false`, it returns `nothing`.
"""
function thinning(chain::Chain)
    @unpack is_ordered, thinning = chain
    is_ordered ? thinning : nothing
end

"""
$(SIGNATURES)

Return the number of samples from the beginning which should be considered warmup.
"""
warmup(chain::Chain) = chain.warmup

"""
$(SIGNATURES)

Return `true` iff the chain is ordered, ie its columns can be considered consecutive draws
from a Markov chain.
"""
is_ordered(chain::Chain) = chain.is_ordered

function Base.vcat(chain1::Chain, chains::Chain...)
    schema1 = schema(chain1)
    @argcheck all(c -> schema(c) == schema1, chains)
    Chain(schema1, mapreduce(sample_matrix, vcat, (chain1, chains...));
          thinning = 1, is_ordered = false, warmup = 0)
end

function Base.getindex(chain::Chain, row_index, name::Symbol)
    getindex(sample_matrix(chain), row_index, getproperty(schema(chain), name))
end

function Base.getindex(chain::Chain, row_index, ::Colon)
    getindex(sample_matrix(chain), row_index, schema(chain))
end

Base.eachcol(chain::Chain) = eachcol(sample_matrix(chain))

end


text:
demo_problems_definitions = [
    (x -> sin(x) + sin(10 // 3 * x), (2.7, 7.5), 5.145735, 6.2173, true),
    (x -> -sum(k * sin((k + 1) * x + k) for k = 1:6), (-10, 10), -6.84128, -1.0991, true),
    (x -> -(16x^2 - 24x + 5) * exp(-x), (1.9, 3.9), 2.868034, :left, false),
    (x -> -x * (1.4 - 3x) * sin(18x), (0, 1.2), 0.969171, 1.1416, true),
    (x -> -(x + sin(x)) * exp(-x^2), (-10, 10), 0.67956, -0.6790, true),
    (x -> sin(x) + sin(10 // 3 * x) + log(x) - 0.84x + 3, (2.7, 7.5), 5.19978, :left, true),
    (x -> -sum(k * cos((k + 1) * x + k) for k = 1:6), (-10, 10), -7.1095, 4.9134, true),
    (x -> sin(x) + sin(2 // 3 * x), (3.1, 20.4), 17.039, :right, true),
    (x -> -x * sin(x), (0, 10), 7.9787, :right, true),
    (x -> 2cos(x) + cos(2x), (-π / 2, 2π), 2.09439, 2π, true), # Handle
    (x -> sin(x)^3 + cos(x)^3, (0, 2π), π, :left, true),
    (x -> -x^(2 // 3) - (1 - x^2)^(1 // 3), (0.001, 0.99), 1 / sqrt(2), :left, false),
    (x -> -exp(-x) * sinpi(2x), (0, 4), 0.224885, 0.7248, true),
    (x -> (x^2 - 5x + 6) / (x^2 + 1), (-5, 5), 2.41422, -0.4145, true),
    (x -> -(x - sin(x)) * exp(-x^2), (-10, 10), 1.195137, -1.1951, true),
    (x -> x * sin(x) + x * cos(2x), (0, 10), 4.79507, 9.2039, true),
    (x -> exp(-3x) - sin(x)^3, (0, 20), 9π / 2, 4.7124, true),
]

"""
    demo_problem(i, thin = false)

Return a demo problem for testing `extrema_enclosure` and
`extrema_series`. Returns `(f, (a, b), fmin, fmax, fabsmin, fabsmax)`
where `f` is the function, `(a, b)` gives the interval to use, `fmin`
and `fmax` are the minimum and maximum value on the interval and
`fabsmin` and `fabsmax` are the minimum and maximum value of the
absolute value of `f`.

If `thin = true` then don't return the full interval `[a, b]` but a
thinner one around the location of the minimum.

The problems are taken from
http://infinity77.net/global_optimization/test_functions_1d.html. The
minima given there are only approximate and not even correct in all
cases. To get high precision values we take an a priori approximation
for the location of the minimum and then we rigorously find a zero of
the derivative close to that approximation.
"""
function demo_problem(i::Integer, thin = false)
    f, (a, b), xmin, xmax, crosses = demo_problems_definitions[i]

    a, b = convert(Arf, a), convert(Arf, b)

    if xmin == :left
        xmin = Arb(a)
    elseif xmin == :right
        xmin = Arb(b)
    else
        xmin = convert(Arb, xmin)
        Arblib.add_error!(xmin, Arb(1e-3))
    end

    if xmax == :left
        xmax = Arb(a)
    elseif xmax == :right
        xmax = Arb(b)
    else
        xmax = convert(Arb, xmax)
        Arblib.add_error!(xmax, Arb(1e-3))
    end

    df(x::Arb) = f(ArbSeries((x, 1), degree = 1))[1]
    df(x::ArbSeries) =
        Arblib.derivative(f(ArbSeries(x.arb_poly, degree = Arblib.degree(x) + 1)))

    if !iszero(Arblib.radref(xmin))
        roots, flags = ArbExtras.isolate_roots(df, Arblib.getinterval(Arf, xmin)...)
        @assert only(flags)
        root_min = ArbExtras.refine_root(df, Arb(only(roots)))
        xmin = Arblib.intersect(Arb((a, b)), root_min)
    end

    if !iszero(Arblib.radref(xmax))
        roots, flags = ArbExtras.isolate_roots(df, Arblib.getinterval(Arf, xmax)...)
        @assert only(flags)
        root_max = ArbExtras.refine_root(df, Arb(only(roots)))
        xmax = Arblib.intersect(Arb((a, b)), root_max)
    end

    if thin
        # Thin interval around roots
        interval_min = Arblib.add_error!(copy(xmin), Arf(1e-2))
        interval_max = Arblib.add_error!(copy(xmax), Arf(1e-2))
        # Intersect with interval
        c_min, d_min = getinterval(Arf, intersect(Arb((a, b)), interval_min))
        c_max, d_max = getinterval(Arf, intersect(Arb((a, b)), interval_max))

        return f, (c_min, d_min), f(xmin)
    else
        absmin = crosses ? zero(Arb) : min(abs(f(xmin)), abs(f(xmax)))
        absmax = max(abs(f(xmin)), abs(f(xmax)))
        return f, (a, b), f(xmin), f(xmax), absmin, absmax
    end
end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03nd(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 20
    Q2 = Array{Float64,1}(undef, NMAX)
    E2 = Array{Float64,1}(undef, NMAX-1)
    E = Array{Float64,1}(undef, NMAX-1)
    Q = Array{Float64,1}(undef, NMAX)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    THETA = parse(Float64, replace(vs[2],'D'=>'E'))
    if ( N<0 || N>NMAX )
        @error "Illegal N=$N"
    end

    vs = String[]
    _isz = N
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    Q[1:_isz] .= parsex.(Float64, vs)

    vs = String[]
    _isz = N-1
    while length(vs) < _isz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    E[1:_isz] .= parsex.(Float64, vs)
    close(f)

    println(io, "J:")
    show(io, "text/plain", Bidiagonal(Q[1:N],E[1:N-1],:U))
    println(io)
    Q2[N] = Q[N]^2
    PIVMIN = Q2[N]
    for i in 1:N-1
        Q2[i] = Q[i]^2
        E2[i] = E[i]^2
        PIVMIN = max( PIVMIN, Q2[i], E2[i] )
    end
    SAFMIN = floatmin(1.0)
    PIVMIN = max( PIVMIN*SAFMIN, SAFMIN )

    NUMSV, INFO = SLICOT.mb03nd!(N, THETA, Q2, E2, PIVMIN)
    @test INFO == 0
    println(io, "NUMSV = $NUMSV")
    println(io, "THETA = $THETA")

end # run_mb03nd()


text:
using ClimateModels, Pkg, Documenter, Test, Suppressor

@testset "cmip" begin
    (mm,gm,meta)=cmip()
    @test isapprox(gm["y"][end],285.71875,atol=1)
end

@testset "ipcc" begin
    (dat, dat1, dat2)=ClimateModels.IPCC_fig1a_read()
    (dat_1b,meta_1b)=ClimateModels.IPCC_fig1b_read()
    (dat2a,dat2b,dat2c)=ClimateModels.IPCC_fig2_read()

    df=IPCC_hexagons()
    clv, ttl, colors=ClimateModels.IPCC_fig3_example(df)

    dat4a=ClimateModels.IPCC_fig4a_read()
    dat4b=ClimateModels.IPCC_fig4b_read()
    dat5=ClimateModels.IPCC_fig5_read()

    @test isa(colors[1],Symbol)
end

tmp=ModelConfig()
show(tmp)
@test isa(tmp,AbstractModelConfig)

tmp=PackageSpec(url="https://github.com/JuliaOcean/AirSeaFluxes.jl")
tmp=ModelConfig(model=tmp)
setup(tmp)
@test clean(tmp)=="no task left in pipeline"

@testset "doctests" begin
    doctest(ClimateModels; manual = false)
end


text:
using PowerModelsACDC
using PowerModels
using Memento
using InfrastructureModels
using JuMP

# Suppress warnings during testing.
Memento.setlevel!(Memento.getlogger(InfrastructureModels), "error")
Memento.setlevel!(Memento.getlogger(PowerModelsACDC), "error")
Memento.setlevel!(Memento.getlogger(PowerModels), "error")

using Ipopt
using SCS
using Cbc
using Juniper

using Test

local_test = false   # as some tests require Mosek, only limited set sent to travis.

ipopt_solver = JuMP.optimizer_with_attributes(Ipopt.Optimizer, "tol" => 1e-6, "print_level" => 0)
scs_solver = JuMP.optimizer_with_attributes(SCS.Optimizer, "verbose" => 0)
cbc = JuMP.optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
juniper = JuMP.optimizer_with_attributes(Juniper.Optimizer, "nl_solver" => ipopt_solver, "mip_solver" => cbc, "time_limit" => 7200)

if local_test == true
    ### ONLY for local testing, not supported intravis due to licensces ##############
    import Gurobi
    import Mosek
    import MosekTools
    # import CPLEX
    gurobi = JuMP.optimizer_with_attributes(Gurobi.Optimizer)
    # cplex = JuMP.optimizer_with_attributes(CPLEX.Optimizer)
    mosek = JuMP.optimizer_with_attributes(Mosek.Optimizer)
    ##############################
end
include("common.jl")
data_dc = build_mn_data("../test/data/tnep/case4_original.m")

@testset "PowerModelsACDC" begin

include("pf.jl")

include("opf.jl")

include("tnep.jl")

end


text:
# Check default
@test flex_data.uncertainty_set.name == :Ellipsoid && flex_data2.uncertainty_set.name == :Ellipsoid
@test FlexibilityAnalysis.getuncertaintyset(m).name == :Ellipsoid

# Check safeguards
@test_throws ErrorException setuncertaintyset(m, :Ellipsoid)
@test_throws ErrorException setuncertaintyset(m, :Ellipsoid, 42)
@test_throws ErrorException setuncertaintyset(m, :PNorm)
@test_throws ErrorException setuncertaintyset(m, :ellipsoid)
@test_throws ErrorException setuncertaintyset(m, :PNorm, 3)
@test_throws ErrorException setuncertaintyset(m, :Hyperbox, [[1; 1];[1; 1]])
@test_throws ErrorException setuncertaintyset(m, :Hyperbox, [[[1; 1]];[[1; 1; 1]]])

# Test each set type
a = 0
covar = [11.11 a a a; a 11.11 a a; a a 11.11 a; a a a 11.11]
setcovariance(m, covar)
@test flex_data.uncertainty_set.name == :Ellipsoid
setuncertaintyset(m2, :PNorm, 1)
@test flex_data2.uncertainty_set.name == :PNorm && flex_data2.uncertainty_set.p == 1
setuncertaintyset(m2, :PNorm, 2)
@test flex_data2.uncertainty_set.name == :PNorm && flex_data2.uncertainty_set.p == 2
setuncertaintyset(m2, :PNorm, Inf)
@test flex_data2.uncertainty_set.name == :PNorm && flex_data2.uncertainty_set.p == Inf
setuncertaintyset(m2, :Hyperbox, [[[1; 1]]; [[1; 1]]])
@test flex_data2.uncertainty_set.name == :Hyperbox && all(flex_data2.uncertainty_set.neg_dev .== ones(2))
setuncertaintyset(m2, :PNorm, 1, only_positive = true)
@test flex_data2.uncertainty_set.name == :PNorm && flex_data2.uncertainty_set.only_positive
setuncertaintyset(m2, :Hyperbox, [[[1; 1]]; [[1; 1]]], only_positive = true)
@test flex_data2.uncertainty_set.name == :Hyperbox && flex_data2.uncertainty_set.only_positive
setuncertaintyset(m2, :Ellipsoid, Matrix(I, flex_data2.numRVs, flex_data2.numRVs), only_positive = true)
@test flex_data2.uncertainty_set.name == :Ellipsoid && flex_data2.uncertainty_set.only_positive


text:
sL = 0.13
sR = -0.27

KitBase.linear(sL, sR)
KitBase.vanleer(sL, sR)
KitBase.vanleer(sL, sR, sR)
KitBase.minmod(sL, sR)
KitBase.minmod(sL, sR, sR)
KitBase.vanalbaba(sL, sR)
KitBase.superbee(0.51, 1.99)
KitBase.superbee(0.49, 2.01)
KitBase.superbee(0.51, 2.01)

wL = rand(3)
wR = randn(3)
dx = 1e-2

KitBase.reconstruct2(sL, sR, dx)
KitBase.reconstruct2(wL, wR, dx)
KitBase.reconstruct2(randn(3, 2), randn(3, 2), dx)
KitBase.reconstruct2(randn(3, 2, 2), randn(3, 2, 2), dx)

KitBase.reconstruct2!(similar(wL), wL, wR, dx)
KitBase.reconstruct2!(randn(3, 2), randn(3, 2), randn(3, 2), dx)
KitBase.reconstruct2!(randn(3, 2, 2), randn(3, 2, 2), randn(3, 2, 2), dx)

wN = rand(3)

KitBase.reconstruct3(sL, rand(), sR, dx, dx, :linear)
KitBase.reconstruct3(sL, rand(), sR, dx, dx, :vanleer)
KitBase.reconstruct3(sL, rand(), sR, dx, dx, :minmod)
KitBase.reconstruct3(wL, wN, wR, dx, dx, :linear)
KitBase.reconstruct3(wL, wN, wR, dx, dx, :vanleer)
KitBase.reconstruct3(wL, wN, wR, dx, dx, :minmod)
KitBase.reconstruct3(randn(3, 2), randn(3, 2), randn(3, 2), dx, dx)
KitBase.reconstruct3(randn(3, 2, 2), randn(3, 2, 2), randn(3, 2, 2), dx, dx)

KitBase.reconstruct3!(similar(wL), wL, wN, wR, dx, dx, :linear)
KitBase.reconstruct3!(similar(wL), wL, wN, wR, dx, dx, :vanleer)
KitBase.reconstruct3!(similar(wL), wL, wN, wR, dx, dx, :minmod)
KitBase.reconstruct3!(randn(3, 2), randn(3, 2), randn(3, 2), randn(3, 2), dx, dx)
KitBase.reconstruct3!(
    randn(3, 2, 2),
    randn(3, 2, 2),
    randn(3, 2, 2),
    randn(3, 2, 2),
    dx,
    dx,
)
KitBase.reconstruct3!(
    randn(3, 2, 2, 2),
    randn(3, 2, 2, 2),
    randn(3, 2, 2, 2),
    randn(3, 2, 2, 2),
    dx,
    dx,
)

KitBase.weno5(-2.0, -1.0, 0.0, 1.0, 2.0)

#--- filter ---#
let deg = 2, u = rand(deg + 1)
    ℓ = rand(deg + 1)

    modal_filter!(u, 1e-6; filter = :l2)
    modal_filter!(u, 1e-6; filter = :l2opt)
    modal_filter!(u, 1e-6, ℓ; filter = :l1)
    modal_filter!(u, ℓ; filter = :lasso)
    modal_filter!(u, 10; filter = :exp)
    modal_filter!(u, 10; filter = :houli)
end

let deg = 2, u = rand(deg + 1, deg + 1)
    modal_filter!(u, 1e-6, 1e-6; filter = :l2)
    modal_filter!(u, 1e-6, 1e-6; filter = :l2opt)
end