text:
###############################################################################
include("copyFolder.jl")
using Base.Filesystem
MAINFolder="---"
try
	global MAINFolder
	MAINFolder = ENV["BACKUPDATAFOLDER"]
	# println( MAINFolder )
catch e
	println("请定义环境变量  BACKUPDATAFOLDER 为   本地主机硬盘存储数据文件的文件夹位置 ，比如： set BACKUPDATAFOLDER=e:\\data")
	exit()
end

println( "Local Machine Folder is:", MAINFolder)

function cpTOLocalDIsk(subfolders= ["b","tem","t"] )
 	srcRoot= "../"
	tgrRoot=MAINFolder
	for folder in subfolders
        proCopyFolder( "$srcRoot/$folder","$tgrRoot/$folder")
	end
end

function cpToUSBDIsk(subfolders= ["b","tem","t"])
	srcRoot=MAINFolder
	tgrRoot= "../"
	for folder in  subfolders
        proCopyFolder( "$srcRoot/$folder","$tgrRoot/$folder")
	end
end

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

subfolders= ["b","tem","t"]
try
	global subfolders
	subfolders1 = ENV["SUBFOLDERS"]
	subfolders=split(subfolders1,":")
	println( "::",subfolders )
catch e
	println("请定义环境变量  SUBFOLDERS 为 从主机复制数据到本盘（toUSBDisk），还是从本盘复制数据到主机（fromUSBDisk） ，比如： set BACKUPDATAFOLDER=e:\\data")
end
cmd_STRING=""
try
	global cmd_STRING
	cmd_STRING=ENV["CMD_STRING"] # toUSBDisk or fromUSBDisk
	# subfolders = ENV["SUBFOLDERS"].split(":")

catch e
	println("请定义环境变量  CMD_STRING 为 从主机复制数据到本盘（toUSBDisk），还是从本盘复制数据到主机（fromUSBDisk） ，比如： set BACKUPDATAFOLDER=e:\\data")
	exit()
end
println("--:", cmd_STRING )
if cmd_STRING=="toUSBDisk"
	cpToUSBDIsk(subfolders)
else
	if cmd_STRING=="fromUSBDisk"
		cpTOLocalDIsk(subfolders)
	else
		println( "::请定义环境变量的值 CMD_STRING 为 fromUSBDisk 或者 toUSBDisk ")
	end
end
#
#


text:
# ---
# title: 1297. Maximum Number of Occurrences of a Substring
# id: problem1297
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Bit Manipulation
# link: <https://leetcode.com/problems/maximum-number-of-occurrences-of-a-substring/description/>
# hidden: true
# ---
#
# Given a string `s`, return the maximum number of ocurrences of **any**
# substring under the following rules:
#
#   * The number of unique characters in the substring must be less than or equal to `maxLetters`.
#   * The substring size must be between `minSize` and `maxSize` inclusive.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aababcaab", maxLetters = 2, minSize = 3, maxSize = 4
#     Output: 2
#     Explanation: Substring "aab" has 2 ocurrences in the original string.
#     It satisfies the conditions, 2 unique letters and size 3 (between minSize and maxSize).
#
#
# **Example 2:**
#
#
#
#     Input: s = "aaaa", maxLetters = 1, minSize = 3, maxSize = 3
#     Output: 2
#     Explanation: Substring "aaa" occur 2 times in the string. It can overlap.
#
#
# **Example 3:**
#
#
#
#     Input: s = "aabcabcab", maxLetters = 2, minSize = 2, maxSize = 3
#     Output: 3
#
#
# **Example 4:**
#
#
#
#     Input: s = "abcde", maxLetters = 2, minSize = 3, maxSize = 3
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 10^5`
#   * `1 <= maxLetters <= 26`
#   * `1 <= minSize <= maxSize <= min(26, s.length)`
#   * `s` only contains lowercase English letters.
#
#
## @lc code=start
using LeetCode

function max_freq(s::String, max_letters::Int, min_size::Int, max_size::Int)::Int
    counts = Dict{SubString,Int}()

    for i in 1:(length(s) - min_size + 1),
        j in (i + min_size - 1):min(length(s), i + max_size + 1)

        ŝ = SubString(s, i, j)
        if length(Set(ŝ)) <= max_letters
            counts[ŝ] = get(counts, ŝ, 0) + 1
        end
    end
    return length(counts) > 0 ? maximum(values(counts)) : 0
end
## @lc code=end

ext: jl


text:
using TimeZones, CSV
using Mise

function read_data()
    input_dir = "/input"
    obs_path = joinpath(input_dir, "station", "aerosol_observatory_2017_aironly.xlsx")
    aes_dir = "aerosol"
    wea_dir = joinpath("weather", "seoul")
    start_date = ZonedDateTime(2008, 1, 1, 1, 0, tz"Asia/Seoul")
    end_date = ZonedDateTime(2020, 10, 31, 23, 0, tz"Asia/Seoul")

    join_data(input_dir, obs_path, aes_dir, wea_dir, start_date, end_date)
end

function read_raw_PM1025()
    input_dir = "/input"
    aes_dir = "aerosol"
    start_date = ZonedDateTime(2008, 1, 1, 1, 0, tz"Asia/Seoul")
    end_date = ZonedDateTime(2020, 10, 31, 23, 0, tz"Asia/Seoul")

    stn_code = 111123
    df_aes = parse_aerosols(aes_dir, input_dir)
    df_aes_jongro = filter_station(df_aes, stn_code)

    # check path
    mkpath(input_dir)
    filename = joinpath(input_dir, "input_raw_PM1025.csv")
    CSV.write(filename, df_aes_jongro)
end

read_data()
#read_raw_PM1025()


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_sastre_opt_m5_rho1_68(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sastre_opt_m5_rho1_68!(A_copy)
end

@inline function exp_sastre_opt_m5_rho1_68!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=8
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 T2k8
    # Computing Bb2 = x*I+x*A
    coeff1=-0.2940691417403946
    coeff2=0.4830197021196582
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=-0.2940691417403946
    coeff2=0.4830197021196582
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=-0.3371652952117645
    coeff2=0.4456299027272243
    coeff3=0.08435039764895878
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.04182430699533517
    coeff2=-0.19033210762698088
    coeff3=0.9365591932640623
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=0.0007801504959769123
    coeff2=0.036934048204913826
    coeff3=-0.33589987345926525
    coeff4=0.9668380398527914
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=0.008399213058466935
    coeff2=0.0017935941807155679
    coeff3=0.0005218334570419188
    coeff4=7.480755260549251e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.11409149134408365
    coeff2=0.68729265604707
    coeff3=0.37829567639304673
    coeff4=0.033903885370555546
    coeff5=0.36476063896683886
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=0.1382844980970843
    coeff2=0.13000362502263652
    coeff3=0.06285653463395623
    coeff4=0.018703010057794192
    coeff5=0.6395504919932988
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.968291126022028
    coeff2=0.5165851751813124
    coeff3=0.47391353450112705
    coeff4=0.11560763975481574
    coeff5=5.160367463472582
    coeff6=0.6818767988800816
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.968291126022028
    coeff2=0.5165851751813124
    coeff3=0.47391353450112705
    coeff4=0.11560763975481574
    coeff5=5.160367463472582
    coeff6=0.6818767988800816
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing T2k8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=-0.055527011448072344
    coeff2=-0.019592818624400106
    coeff3=-0.01516878804519063
    coeff4=-0.002839711497134736
    coeff5=0.014357966768578765
    coeff6=0.009741431366481671
    coeff7=1.0239727668224385
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    return memslots1 # Returning T2k8
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
include("../aimajulia.jl");

using Test;

using Main.aimajulia;

#The following mdp tests are from the aima-python doctests

tm = Dict([Pair("A", Dict([Pair("a1", (0.3, "B")), Pair("a2", (0.7, "C"))])),
            Pair("B", Dict([Pair("a1", (0.5, "B")), Pair("a2", (0.5, "A"))])),
            Pair("C", Dict([Pair("a1", (0.9, "A")), Pair("a2", (0.1, "B"))]))]);

mdp = MarkovDecisionProcess("A", Set(["a1", "a2"]), Set(["C"]), tm, states=Set(["A","B","C"]));

@test (transition_model(mdp, "A", "a1") == (0.3, "B"));

@test (transition_model(mdp, "B", "a2") == (0.5, "A"));

@test (transition_model(mdp, "C", "a1") == (0.9, "A"));

@test (repr(value_iteration(aimajulia.sequential_decision_environment, epsilon=0.01)) ==
        "Dict((2, 3)=>0.486437,(2, 1)=>0.398102,(3, 1)=>0.509285,(1, 4)=>0.129589,(3, 3)=>0.795361,(1, 3)=>0.344613,(3, 2)=>0.649581,(2, 4)=>-1.0,(1, 1)=>0.295435,(1, 2)=>0.253487,(3, 4)=>1.0)");

pi = optimal_policy(aimajulia.sequential_decision_environment, value_iteration(aimajulia.sequential_decision_environment, epsilon=0.01));

@test (repr(pi) == "Dict{Any,Any}((2, 3)=>(1, 0),(2, 1)=>(1, 0),(3, 1)=>(0, 1),(1, 4)=>(0, -1),(3, 3)=>(0, 1),(1, 3)=>(1, 0),(3, 2)=>(0, 1),(2, 4)=>nothing,(1, 1)=>(1, 0),(1, 2)=>(0, 1),(3, 4)=>nothing)");

@test (repr(to_arrows(aimajulia.sequential_decision_environment, pi)) == "Union{Nothing, String}[\"v\" \">\" \"v\" \"<\"; \"v\" nothing \"v\" \".\"; \">\" \">\" \">\" \".\"]");

@test (policy_iteration(aimajulia.sequential_decision_environment) ==
        Dict([Pair((2,3),(1,0)),
            Pair((2,1),(1,0)),
            Pair((3,1),(0,1)),
            Pair((1,4),(0,-1)),
            Pair((3,3),(0,1)),
            Pair((1,3),(1,0)),
            Pair((3,2),(0,1)),
            Pair((2,4),nothing),
            Pair((1,1),(1,0)),
            Pair((1,2),(0,1)),
            Pair((3,4),nothing)]));

max_stars_repo_name: mikhail-j/aima-julia-aimacode-fork

max_issues_repo_name: mikhail-j/aima-julia-aimacode-fork

max_forks_repo_name: mikhail-j/aima-julia-aimacode-fork


text:
module opentelemetry
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
module proto
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
module common
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
module v1
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
include("common_pb.jl")
end
end
module resource
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
module v1
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
include("resource_pb.jl")
end
end
module trace
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
module v1
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
include("trace_pb.jl")
end
end
module logs
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
module v1
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
include("logs_pb.jl")
end
end
module metrics
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
module v1
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
include("metrics_pb.jl")
end
end
module collector
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
module trace
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
module v1
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
include("trace_service_pb.jl")
end
end
module metrics
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
module v1
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
include("metrics_service_pb.jl")
end
end
module logs
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
module v1
const _ProtoBuf_Top_ = @static if isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_)
    (parentmodule(@__MODULE__))._ProtoBuf_Top_
else
    parentmodule(@__MODULE__)
end
include("logs_service_pb.jl")
end
end
end
end
end


text:
input = readlines("./input")
nums = sort(parse.(Int, input))
set = Set(nums)

for i in 1:length(nums)
  X = nums[i]
  for j in (i + 1):length(nums)
    Y = nums[j]
    target = 2020 - X - Y
    if target in set
      println(X * Y * target)
      exit(0)
    end
  end
end


text:
"""
For first metadata pass we need already cut the area o boolean arrays that we are intrested in and portion of metadata array that we are intrested in
this way indexing will be simpler we will start from 0 and we will reduce memor usage

3) first metadata pass we add to the metadata offset where each block will put its main results and padding results - so all will be stored in result quueue but in diffrent spots
    we need to make ques for paddings longer than number of possible results becouse of possible modifications from neighbouring blocks that can happen simultanously
    we will establish this offsets using atomics- at this pass we will also prepare first work queue with indicies of metadata blocks and booleans indicating is it related to gold pass dilatation step or other pass

5) we do the metadata pass we analyze only those blocks that are in the borders of intrests - max min x y z was specified in step 1 , we check weather block is set to be activated
 and in case it is not full we will make it active , if the block is set as active we just add it to work queue  that is appropriate to next iteration
  - we need also to scan  the border result ques if there is any duplicate result - if so we set it to zeros and we reduce the border result counter
        of course we check it only in case new counter is bigger than old counter

        so 1) we check metadata when metadata comply with our predicate described above we atomically increase local shared memory work queue counter
we synchronize
        -we proceed if local workqueue counter is greater than 0
        2)  we add the local workqueue counter to global keep old value as an offset in shared memory
        3) sync threads,  add to the work queue the data from registers of those threads that met the predicate
"""

module  MetadataAnalyzePass
using CUDA, Logging,..CUDAGpuUtils, ..ResultListUtils,..WorkQueueUtils,..ScanForDuplicates, Logging,StaticArrays, ..IterationUtils, ..ReductionUtils, ..CUDAAtomicUtils,..MetaDataUtils, ..ScanForDuplicates
export @metaDataWarpIter, @loadCounters,@analyzeMetadataFirstPass, @checkIsActiveOrFullOr,@setIsToBeActive

"""
this will enable iteration of metadata block
linear index is the same for each threadIdxX in a block hence the bigger in y direction is thread block the more threads will work on single metadata block

"""
macro metaDataWarpIter(metaDataDims,loopWarpMeta,metaDataLength,ex)

    return  esc(quote
    isInRange = true
    linIdexMeta = UInt32(0)

    # gridDimX()*blockDimX()

    @unroll for j in 0:($loopWarpMeta-1)
        linIdexMeta= ((blockIdxX()-1)*$loopWarpMeta*blockDimX())+j*blockDimX()+ threadIdxX()-1
        # linIdexMeta= ((blockIdxX()-1)*(gridDimX())*blockDimX()*$loopWarpMeta)+j*blockDimX()+ threadIdxX()
        isInRange=(linIdexMeta<$metaDataLength)
        xMeta= rem(linIdexMeta,$metaDataDims[1])
        zMeta= fld((linIdexMeta),$metaDataDims[1]*$metaDataDims[2])
        yMeta= fld((linIdexMeta-((zMeta*$metaDataDims[1]*$metaDataDims[2] ) + xMeta )),$metaDataDims[1])
        # zMeta= fld(zMeta,32)
      $ex
    end

    end)

    # mainExp = generalizedItermultiDim(;xname=:(xMeta)
    # ,yname= :(yzSpot)
    # ,arrDims=metaDataDims
    # ,loopXdim=loopXMeta
    # ,loopYdim=loopYZMeta
    #  ,yOffset = :(ydim*gridDim().x)
    # ,yAdd=  :(blockIdxX()-1)
    # ,additionalActionBeforeY= :( yMeta= rem(yzSpot,$metaDataDims[2]) ; zMeta= fld(yzSpot,$metaDataDims[2]) )
    # ,additionalActionBeforeX= :( isInRange = ( yMeta < $metaDataDims[2] && zMeta<$metaDataDims[3] && xMeta <= $metaDataDims[1]  ) )
    # ,nobundaryCheckX=true
    # , nobundaryCheckY=true
    # , nobundaryCheckZ =true
    # # ,yCheck = :(yMeta < $metaDataDims[2] && zMeta<$metaDataDims[3] )
    # # ,xCheck = :(xMeta <= $metaDataDims[1])
    # # ,xAdd= :(threadIdxX()-1)# to keep all 0 based
    # ,is3d = false
    # , ex = ex)
    # return esc(:( $mainExp))
end

"""
now we upload all data related to amount of data that is of our intrest
as we need to perform basically the same work across all warps - instead on specializing threads in warp we will execute the same fynction across all warps
so warp will execute the same function just with varying data as it should be

"""
macro loadCounters()
    return esc(quote
        @unroll for i in 1:14
           @exOnWarp i begin
            if(isInRange)
                shmemSum[threadIdxX(),i]= @accMeta(getBeginingOfFpFNcounts()+ i)
            end
           end
        end#for

        # now in order to get offsets we need to atomially access the resOffsetCounter - we add to them total fp or fn cout so next blocks will not overwrite the
        # area that is scheduled for this particular metadata block
        # we need to supply linear coordinate for atomicallyAddToSpot
        @exOnWarp 15 begin
            if(isInRange)
            count = @accMeta(getBeginingOfFpFNcounts()+ 15)
                if(count>0)
                    shmemPaddings[threadIdxX()]= true
                    shmemSum[threadIdxX(),15]= atomicAdd(globalFpResOffsetCounter,  ceil(count*1.5)  )+1
                else
                    shmemSum[threadIdxX(),15]= 0
                end
            end
        end

        @exOnWarp 16 begin
            if(isInRange)
            count = @accMeta(getBeginingOfFpFNcounts()+ 16)
                if(count>0)
                    shmemSum[threadIdxX(),16]= atomicAdd(globalFnResOffsetCounter,  ceil( count*1.5 )  )+1
                    shmemPaddings[threadIdxX()+32]= true
                else
                    shmemSum[threadIdxX(),16]= 0
                end
            end
        end

        end)#quote
end #loadCounters

 """
 we analyze metadate as described above
 minX, minY,minZ - minimal indexes of metadata that holds all of the data that is of intrest to us
 maxX,maxY,maxZ  - maximal indexes of metadata that holds all of the data that is of intrest to us
 metaData - global memory data structure that we analyze
 shmemSum - shared memory used primary for reductions
 globalFpResOffsetCounter, globalFnResOffsetCounter  - counters accessed atomically that points where we want to set the  results from this metadata block
 workQueaueA, workQueaueAcounter
 tobeEx - some register boolean that we will reuse
 """
 macro analyzeMetadataFirstPass()
         return esc(quote
         #clear shmem paddings
        @ifY 1 shmemPaddings[threadIdxX()]=false
        @ifY 2 shmemPaddings[threadIdxX()+32]=false
        @ifXY 3 3 workCounterLocalInShmem[1]=0
        sync_threads()
         # we need to iterate over all metadata blocks with checks so the blocks can not be  outside the area of intrest defined by  minX, minY,minZ and maxX,maxY,maxZ
         @metaDataWarpIter(metaDataDims,loopWarpMeta,metaDataLength,begin
             #now we upload all data related to amount of data that is of our intrest
             #as we need to perform basically the same work across all warps - instead on specializing threads in warp we will execute the same fynction across all warps
             # so warp will execute the same function just with varying data as it should be
            # if(isInRange)
            #   CUDA.@cuprint "xMeta $(xMeta) yMeta $(yMeta) zMeta $(zMeta) isInRange $(isInRange) linIdexMeta $(linIdexMeta) metaDataDims $(metaDataDims) loopWarpMeta $(loopWarpMeta) metaDataLength $(metaDataLength) \n "
            # end
            @loadCounters()

            sync_threads()

             ######  we need to establish is block is active at the first pass block is active simply  when total count of fp and fn is greater than 0
            #we are adding 1 to meta y z becouse those are 0 based ...

            @ifY 1 if(shmemPaddings[threadIdxX()+32] && isInRange) begin
                @appendToLocalWorkQueue((xMeta+1),(yMeta+1),(zMeta+1),0 )
                end
            end
            @ifY 2 if(shmemPaddings[threadIdxX()] && isInRange) begin
                 @appendToLocalWorkQueue((xMeta+1),(yMeta+1),(zMeta+1),1 )
                    end
                end

            @exOnWarp 3 if((shmemSum[threadIdxX(),15]) >0 && isInRange) setBlockasCurrentlyActiveInSegm(metaData, xMeta+1,yMeta+1,zMeta+1)    end
            @exOnWarp 4 if((shmemSum[threadIdxX(),16]) >0 && isInRange) setBlockasCurrentlyActiveInGold(metaData, xMeta+1,yMeta+1,zMeta+1)     end

            #####3set offsets
             #now we will calculate and set the result queue offsets for each offset we need to synchronize warps in order to have unique offsets
             #we can not parallalize it more as we need to sequentially set offsets

            @exOnWarp 5 begin if((shmemSum[threadIdxX(),15]) >0 && isInRange) @unroll for i in 0:6
                     #set fp
                   value=floor(shmemSum[threadIdxX(),15])+1
                   if(i>0)
                    value+= ceil(shmemSum[threadIdxX(),((i-1)*2+1)]*1.45)
                   end
                   shmemSum[threadIdxX(),15]= value
                   @setMeta(((getResOffsetsBeg()-1) +i*2+1) ,value)
                     end#for
                    end #if
                end

            @exOnWarp 6 begin if((shmemSum[threadIdxX(),16]) >0 && isInRange) @unroll for i in 0:6
                 #set fn
                 value=shmemSum[threadIdxX(),16]
                 if(i>0)
                    value+= ceil(shmemSum[threadIdxX(),((i-1)*2+2)]*1.45)+1 #multiply as we can have some entries potentially repeating
                 end
                 shmemSum[threadIdxX(),16]= value
                 @setMeta(((getResOffsetsBeg()-1) +i*2+2) ,value)
                end#for
            end#if
        end

        #so now we need to make sure that we have enough space left in the shmemblockData for work queue elements
        sync_threads()
        if(workCounterLocalInShmem[1]+32>= shmemblockDataLenght)
            @appendToGlobalWorkQueue()
            sync_threads()
            @ifXY 1 1 workCounterLocalInShmem[1]=0
            sync_threads()
        end

         end)# outer loop expession  )
         # probably we do not need to clear as we assign not adding values ...
         #clearSharedMemWarpLong(shmemSum, UInt8(14), Float32(0.0))

         sync_threads()
         @appendToGlobalWorkQueue()
        #  @ifXY 1 1 CUDA.@cuprint "  aaaa workQueueCounter[1] $(workQueueCounter[1]) \n"

        end )
 end
"""
establish is the  block  is active full or be activated, and we are saving this information into surcehmem
"""
macro checkIsActiveOrFullOr()
    return esc(quote
        @exOnWarp 30 if(isInRange) shmemPaddings[(threadIdxX())] = @accMeta(getFullInGoldNumb() ) end#  isBlockFulliInGold(metaData, xMeta,yMeta+1,zMeta+1)
        @exOnWarp 31 if(isInRange)   shmemPaddings[(threadIdxX())+33] = @accMeta(getIsToBeActivatedInGoldNumb() ) end # isBlockToBeActivatediInGold(metaData, xMeta,yMeta+1,zMeta+1)
        @exOnWarp 32 if(isInRange)  shmemPaddings[(threadIdxX())+33*2] = @accMeta(getActiveGoldNumb() ) end # isBlockCurrentlyActiveiInGold(metaData, xMeta,yMeta+1,zMeta+1)

        @exOnWarp 33 if(isInRange) shmemPaddings[(threadIdxX())+33*3] =@accMeta(getFullInSegmNumb()) end # isBlockFullInSegm(metaData, xMeta,yMeta+1,zMeta+1)
        @exOnWarp 34 if(isInRange) shmemPaddings[(threadIdxX())+33*4] = @accMeta(getIsToBeActivatedInSegmNumb() ) end # isBlockToBeActivatedInSegm(metaData, xMeta,yMeta+1,zMeta+1)
        @exOnWarp 35 if(isInRange) shmemPaddings[(threadIdxX())+33*5] = @accMeta(getActiveSegmNumb()) end # isBlockCurrentlyActiveInSegm(metaData, xMeta,yMeta+1,zMeta+1)
end)#quote
end#checkIsActiveOrFullOr

"""
given data in sourceShmem loaded by checkIsActiveOrFullOr() we will  mark the block as active  ( or not)
    and if is to be active add it to work queue
"""
macro setIsToBeActive()
    return esc(quote
    # CUDA.@cuprint " a $(!(shmemPaddings[(threadIdxX())])) b $(((shmemPaddings[(threadIdxX())+33]))) c $(shmemPaddings[(threadIdxX())+33*2])) \n"

    #if(isInRange) CUDA.@cuprint " a $(!(shmemPaddings[(threadIdxX())])) b $(((shmemPaddings[(threadIdxX())+33]))) c $(shmemPaddings[(threadIdxX())+33*2])) \n" end

        @exOnWarp 1 if((!(shmemPaddings[(threadIdxX())]) && ((shmemPaddings[(threadIdxX())+33]) ||  (shmemPaddings[(threadIdxX())+33*2]))) &&isInRange  )
                # CUDA.@cuprint "set meta other pass \n  "
                @setMeta(getActiveGoldNumb(),1)
                        @appendToLocalWorkQueue(xMeta+1,yMeta+1,zMeta+1,1 )
                    end
        @exOnWarp 2 if((!(shmemPaddings[(threadIdxX())+33*3])  && ((shmemPaddings[(threadIdxX())+33*4])  ||  (shmemPaddings[(threadIdxX())+33*5]))) &&isInRange )
                    # CUDA.@cuprint "set meta other pass \n  "

                    @setMeta(getActiveSegmNumb(),1)
                        @appendToLocalWorkQueue(xMeta+1,yMeta+1,zMeta+1,0 )
                    end
    end)#quote

end

    """
    will be invoked in order to iterate over the metadata  after some dilatations were already done - we need to
        establish is block to be activated or inactivated or left as is
        if block is active it needs to be added to work queue
        using some spare threads we will also housekeeping like for example switching active work queue etc
        we will check rescounters of border res ques and compare with old ones - if any will be grater than old we will scan for any repeating results
            - it could be the case that neighbouring blocks concurently added the same results - in this case we need to set one of those to 0 and reduce the counter
        we will do all by using single warp per metadata block
        globalCurrentFpCount, globalCurrentFnCount - representing current number of already covere fp and fns
    """
    macro setMEtaDataOtherPasses(locArr,offsetIter,iterThrougWarNumb)
        return esc(quote
        #clear
        $locArr=0
        $offsetIter=0
        isMaskFull=false
    for iterY in 0:inBlockLoopXZIterWithPadding
        if((threadIdxY()+iterY*dataBdim[2])<=dataBdim[1]  )
            #we are reusing offsetIter
            for bitPos in 1:7
                shmemPaddings[threadIdxX(),(threadIdxY()+iterY*dataBdim[2]),bitPos ]=false
            end
        end
    end
    sync_threads()

        @metaDataWarpIter(metaDataDims,loopWarpMeta,metaDataLength,begin
        isMaskOkForProcessing=false
            #first we will check is block full active or be activated and we will set later on this basis what blocks should be put to work queue
             @checkIsActiveOrFullOr()

            #now we need to go through  those numbers and in case some of the border queues were incremented we need to analyze those added entries to establish is there
            # any duplicate in case there will be we need to decrement counter and set the corresponding duplicated entry to 0
            #here we load data about wheather there is anything to be validated here - we save data so it can be read from the perspective of this block
            # #and the blocks aroud that will want to analyze paddings
            @loadAndScanForDuplicates(iterThrougWarNumb,$locArr,$offsetIter)

        # #     #we set information that block should be activated in gold  and segm
             @setIsToBeActive()
            #just to be sure that we will have enough space for workque items
             sync_threads()
             if(workCounterLocalInShmem[1]+32>= shmemblockDataLenght)
                 @appendToGlobalWorkQueue()
                 sync_threads()
                 @ifXY 1 1 workCounterLocalInShmem[1]=0
                 sync_threads()
             end

        end    )
        sync_threads()

        #now we add to the global variables all of the fps and fns after corrections for duplicates
        @ifXY 1 1 begin
            atomicAdd(globalCurrentFpCount, alreadyCoveredInQueues[1]+ alreadyCoveredInQueues[3]+ alreadyCoveredInQueues[5]+ alreadyCoveredInQueues[7]+ alreadyCoveredInQueues[9]+ alreadyCoveredInQueues[11]+ alreadyCoveredInQueues[13])
        end
            @ifXY 2 1 atomicAdd(globalCurrentFnCount, alreadyCoveredInQueues[2]+ alreadyCoveredInQueues[4]+ alreadyCoveredInQueues[6]+ alreadyCoveredInQueues[8]+ alreadyCoveredInQueues[10]+ alreadyCoveredInQueues[12]+ alreadyCoveredInQueues[14])

            sync_threads()

            @metaDataWarpIter(metaDataDims,loopWarpMeta,metaDataLength,begin
            #now we need to set old caounters to the value of new counters so at next dilatation we will count only new values ...
            for i in 1:14
                if(isInRange)
                  @exOnWarp (i+37) @setMeta((getOldCountersBeg() +i),@accMeta(getNewCountersBeg() +i))
                end
            end
            end)
            #  ########clear
             for i in 1:14
                @exOnWarp (i+23) shmemSum[threadIdxX(),i]= 0
             end
            # @iterateLinearly shmemblockDataLoop shmemblockDataLenght shmemblockData[i]=0
            @appendToGlobalWorkQueue()
            $locArr=0
            $offsetIter=0

         end )
    end

  end #MetadataAnalyzePass


text:
using Plots
using OceanParameterizations

files =  ["free_convection", "strong_wind", "strong_wind_no_coriolis", "weak_wind_strong_cooling",
          "strong_wind_weak_cooling", "strong_wind_weak_heating"]

## uw profiles unscaled

Ts3 = Dict() # maps file name to T array
for file in files
    𝒟 = data(file, reconstruct_fluxes=true) # <: OceananigansData
    Ts3[file] = 𝒟.uw.coarse
end

Ts3

anim = @animate for i in 1:288
    plot(legend=false, xlabel="uw", ylabel="z", xlims=(-0.0006,0.0001))
    for (file, T) in Ts3
        plot!(T[:,i], collect(1:33), title = "Unscaled", label="$(file)")
    end
end
gif(anim, pwd()*"/uw_reconstructed_with_subgrid.gif", fps=20)

## vw profiles unscaled

Ts3 = Dict() # maps file name to T array
for file in files
    𝒟 = data(file, reconstruct_fluxes=true) # <: OceananigansData
    Ts3[file] = 𝒟.vw.coarse
end

Ts3

anim = @animate for i in 1:288
    plot(legend=false, xlabel="vw", ylabel="z", xlims=(-0.0001,0.0005))
    for (file, T) in Ts3
        plot!(T[:,i], collect(1:33), title = "Unscaled", label="$(file)")
    end
end
gif(anim, pwd()*"/vw_reconstructed_with_subgrid.gif", fps=20)

## wT profiles unscaled

Ts3 = Dict() # maps file name to T array
for file in files
    𝒟 = data(file, reconstruct_fluxes=true) # <: OceananigansData
    Ts3[file] = 𝒟.wT.coarse
end

Ts3

anim = @animate for i in 1:288
    plot(legend=false, xlabel="wT", ylabel="z", xlims=(-0.00009,0.00006))
    for (file, T) in Ts3
        plot!(T[:,i], collect(1:33), title = "Unscaled", label="$(file)")
    end
end
gif(anim, pwd()*"/wT_reconstructed_with_subgrid.gif", fps=20)


text:
using Documenter, POMDPModelTools

makedocs(
    modules = [POMDPModelTools],
    format = :html,
    sitename = "POMDPModelTools.jl"
)

deploydocs(
    repo = "github.com/JuliaPOMDP/POMDPModelTools.jl.git",
    julia = "1.0",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
#! /usr/bin/env julia

using BoundingBoxes
using Base.Test

# Initialize Bounding box
@boundingbox Bounds2F, "x", "y"
bb2f = Bounds2F{Float64}(1,2,3,4)

@test bb2f.x_max == 1
@test bb2f.y_max == 2
@test bb2f.x_min == 3
@test bb2f.y_min == 4

@boundingbox Bounds2I, "x", "y"
bb2i = Bounds2I{Int64}(4,3,2,1)

@test bb2i.x_max == 4
@test bb2i.y_max == 3
@test bb2i.x_min == 2
@test bb2i.y_min == 1

# test equality
bb2e = Bounds2I{Int64}(4,3,2,1)

@test bb2i == bb2e
@test bb2i != bb2f

# test empty bounds
bb2f = Bounds2I(Float64)

@test bb2f == Bounds2I{Float64}(-Inf,-Inf,Inf,Inf)
@test typeof(bb2f) == Bounds2I{Float64}

# test arg promoting
bb2a = Bounds2I(1,2,3,4.0)
@test typeof(bb2a) == Bounds2I{Float64}

bb2b = Bounds2I(1,2,3,4.0+im)
@test typeof(bb2b) == Bounds2I{Complex{Float64}}

# test bounds updating
@boundingbox Bounds3, "x", "y", "z"
bb3a = Bounds3(Float64)
@test update!(bb3a, [1,2,3]) == nothing
update!(bb3a, [4,5,6])
@test bb3a.x_max == 4
@test bb3a.y_max == 5
@test bb3a.z_max == 6
@test bb3a.x_min == 1
@test bb3a.y_min == 2
@test bb3a.z_min == 3

# test in
bb3a = Bounds3{Float64}(0.5,0.5,0.5,0.25,0.25,0.25)
bb3b = Bounds3{Float64}(1,1,1,0,0,0)
@test in(bb3a,bb3b)
@test !in(bb3b,bb3a)


text:
using HDF5, MPI, Dates
using CodedComputing

function write_coded_input(X; inputfile::AbstractString, name="C", coderate=1)
    nc = ceil(Int, size(X, 1) / coderate) # number of coded rows
    C = Matrix(encode_hadamard(X, nc))
    h5write(inputfile, name, C)
    return C
end

function benchmark_main(directory=joinpath("./simulations/1000genomes", "$(now())"))

    # MNIST
    # directory = "/home/albin/.julia/dev/CodedComputing/simulations/mnist/2021-01-12T10:34:57.725"

    # 1000 genomes
    # directory = "/home/albin/.julia/dev/CodedComputing/simulations/1000genomes/2021-01-12T17:16:23.558"
    directory = "/home/albin/.julia/dev/CodedComputing/simulations/1000genomes/2021-01-14T13:07:35.092"

    # prepare a directory for storing everything
    mkpath(directory)
    inputfile = joinpath(directory, "inputfile.h5")

    # simulation parameters
    kernel = "./src/pca/pca.jl"
    inputdataset = "X"
    nworkers = 42
    ncomponents = 3
    niterations = 5

    iteratedataset = "V0"

    # don't do anything if the input file already exists
    if isfile(inputfile)
        @assert HDF5.ishdf5(inputfile)
    else
        # MNIST dataset
        # X = load_mnist_training()
        # h5write(inputfile, "X", X)

        # 1000 genomes dataset
        X = load_genome_data()
        h5writecsc(inputfile, "X", X)

        V0 = randn(size(X, 2), ncomponents)
        orthogonal!(V0)
        h5write(inputfile, "V0", V0)
    end

    # Measure iteration time
    stepsize = 1.0
    # nwait = nworkers - 2
    nsubpartitions = 1

    for nreplicas in [2, 3, 1]
        npartitions = div(nworkers, nreplicas) * nsubpartitions
        for nwait in [1, npartitions-2, npartitions-1, npartitions] # 3
            # for nwait in [nworkers-10, nworkers-2, nworkers-1, nworkers]
            # for stepsize in [0.5, 0.9, 1.0]

            # with variance reduction
            outputfile = joinpath(directory, "output_$(now()).h5")
            mpiexec(cmd -> run(```
                $cmd -n $(nworkers+1) julia --project $kernel
                $inputfile $outputfile
                --iteratedataset $iteratedataset
                --niterations $niterations
                --nwait $nwait
                --nreplicas $nreplicas
                --ncomponents $ncomponents
                --stepsize $stepsize
                --nsubpartitions $nsubpartitions
                --variancereduced
                --saveiterates
                ```))

            # without variance reduction (regular SGD)
            outputfile = joinpath(directory, "output_$(now()).h5")
            mpiexec(cmd -> run(```
                $cmd -n $(nworkers+1) julia --project $kernel
                $inputfile $outputfile
                --iteratedataset $iteratedataset
                --niterations $niterations
                --nwait $nwait
                --nreplicas $nreplicas
                --ncomponents $ncomponents
                --stepsize $stepsize
                --nsubpartitions $nsubpartitions
                --saveiterates
                ```))
        end
    end

    df = aggregate_benchmark_data(dir=directory, inputfile=inputfile, inputname="X", prefix="output", dfname="df.csv")
end


text:
using SafeTestsets

const GROUP = get(ENV, "GROUP", "All")

if GROUP == "All" || GROUP == "Interface"
  @time @safetestset "AD Tests" begin include("interface/ad.jl") end
  @time @safetestset "Backwards Tests" begin include("interface/backwards.jl") end
  @time @safetestset "Composite Solution Tests" begin include("interface/composite_solution.jl") end
  @time @safetestset "Constrained Time Steps Tests" begin include("interface/constrained.jl") end
  @time @safetestset "Dependent Delay Tests" begin include("interface/dependent_delays.jl") end
  @time @safetestset "Discontinuity Tests" begin include("interface/discontinuities.jl") end
  @time @safetestset "Export Tests" begin include("interface/export.jl") end
  @time @safetestset "Fixed-point Iteration Tests" begin include("interface/fpsolve.jl") end
  @time @safetestset "History Function Tests" begin include("interface/history_function.jl") end
  @time @safetestset "Jacobian Tests" begin include("interface/jacobian.jl") end
  @time @safetestset "Parameterized Function Tests" begin include("interface/parameters.jl") end
  @time @safetestset "saveat Tests" begin include("interface/saveat.jl") end
  @time @safetestset "save_idxs Tests" begin include("interface/save_idxs.jl") end
  @time @safetestset "Unconstrained Time Steps Tests" begin include("interface/unconstrained.jl") end
  @time @safetestset "Units Tests" begin include("interface/units.jl") end
end

if GROUP == "All" || GROUP == "Integrators"
  @time @safetestset "Cache Tests" begin include("integrators/cache.jl") end
  @time @safetestset "Event Tests" begin include("integrators/events.jl") end
  @time @safetestset "Iterator Tests" begin include("integrators/iterator.jl") end
  @time @safetestset "Reinitialization Tests" begin include("integrators/reinit.jl") end
  @time @safetestset "Residual Control Tests" begin include("integrators/residual_control.jl") end
  @time @safetestset "Return Code Tests" begin include("integrators/retcode.jl") end
  @time @safetestset "Rosenbrock Tests" begin include("integrators/rosenbrock.jl") end
  @time @safetestset "SDIRK Tests" begin include("integrators/sdirk.jl") end
  @time @safetestset "Unique Timepoints Tests" begin include("integrators/unique_times.jl") end
  @time @safetestset "Verner Tests" begin include("integrators/verner.jl") end
end

if GROUP == "All" || GROUP == "Regression"
  @time @safetestset "Inference Tests" begin include("regression/inference.jl") end
  @time @safetestset "Waltman Problem Tests" begin include("regression/waltman.jl") end
end


text:
using BenchmarkTools
using JLD

res = JLD.load("result_SmallBenchmark.jld")

for (t,d,l) in [("#### MSA module","MIToS.MSA.","small_bench_msa"),
              ("#### PDB module","MIToS.PDB.","small_bench_mitos_pdb"),
              ("#### Information module","MIToS.Information.","small_bench_information")]
    println("  ")
    println(t, "  ")
    for (k,v) in res["small_result"][l]
        println("  ")
        println("| ", k, " | |  ")
        println("|---|---|  ")
        for (kk,vv) in minimum(v)
            println("| ", replace(kk,d => "") , " | ", BenchmarkTools.prettytime(time(vv)), " |  ")
        end
    end
end

println("#### Pfam/pipeline  ")
println("  ")
println("|   |   |  ")
println("|---|---|  ")
for (k,v) in minimum(res["small_result"]["small_bench_pfam"])
    println("| ", replace(k,"MIToS.Pfam." => "") , " | ", BenchmarkTools.prettytime(time(v)), " |  ")
end

max_stars_repo_name: diegozea/mitos-benchmarks

max_issues_repo_name: diegozea/mitos-benchmarks

max_forks_repo_name: diegozea/mitos-benchmarks


text:
using Test, SafeTestsets

@time @safetestset "planar gradient" begin include("planar.jl") end
@time @safetestset "rand" begin include("rand.jl") end


text:
module ExtraMessages2
include("apriltag_ros.jl")
export apriltag_ros
end


text:
using Dates
using ArgParse
using Printf
using Random
using GaussianProcesses

# before_all = Dates.now()

include("../../covariatefun.jl")
include("../../btg.jl")
path = "../../datasets/"
include(path * "load_unrate.jl")
include(path * "load_unrate_full.jl")

data, label = unrate_full(path)

#data, label = unrate(path)
data, label = unrate_full(path)
max_label = maximum(label)
label = label ./ max_label

xs = copy(data)
xs = xs/maximum(data)

###
###  Manually enter parsed args entries so we can run this file as a script
###
##
parsed_args = Dict()
push!(parsed_args, "randseed" => 4321)
push!(parsed_args, "lmax" => 0.3)
push!(parsed_args, "lmin" => 0.2)
push!(parsed_args, "percent_train" => 0.4)
push!(parsed_args, "p" => 1)

randseed = parsed_args["randseed"]; rng = MersenneTwister(randseed)
Random.seed!(randseed);

# training set
p_train = parsed_args["percent_train"]
n_train = (Int64)(round(p_train * length(data)))
n_test = length(data) - n_train
"""
Input percent train
"""
function sample(inds, percent)
    num = (Int64)(round(percent * length(data)))
    #train = sort(rand(1:490, (num,)))
    ind_train = sort(randperm(length(data))[1:num])
    ind_test = copy(inds)
    deleteat!(ind_test, ind_train)
    #work with indices up until this point
    x_train = xs[ind_train]
    x_test = xs[ind_test]
    return ind_train, x_train, ind_test, x_test
end

ind_train, x_train, ind_test, x_test  = sample(data, p_train) #training and testing pts
train_label = label[ind_train]
test_label = label[ind_test]
y_train = train_label

x = reshape(x_train, n_train, 1)
Fx = covariate_fun(x, parsed_args["p"])

trainingData0 = trainingData(x, Fx, train_label)

d = getDimension(trainingData0); n = getNumPts(trainingData0); p = getCovDimension(trainingData0)

x_test = reshape(x_test, n_test, 1)
Fx_test = covariate_fun(x_test, parsed_args["p"])
y_test = test_label

elapsedmin = 0

#parameter setting
# myquadtype = parsed_args["sparse"] ? ["SparseCarlo", "SparseCarlo"] : ["QuasiMonteCarlo", "QuasiMonteCarlo"]
myquadtype = ["Gaussian", "Gaussian"]
#rangeλ = [-3.0 3]
#rangeλ = reshape([-0.2], 1, 1)
#rangeλ = reshape([-0.2 0.2], 1, 2) #compare to this
#rangeλ = reshape([-0.2], 1, 1) #BETTER BY 0.5 PERCENT
#rangeλ = reshape([-0.001], 1, 1) #compare above to this
rangeλ = reshape([0.1], 1, 1) #compare above to this

#rangeλ = reshape([-2.0 2.0], 1, 2) #BETTER BY 0.5 PERCENT
lmin = parsed_args["lmin"]
lmax = parsed_args["lmax"]
#rangeθ = [1/lmax^2 1/lmin^2]
#rangeθ = reshape([50.0 150.0], 1, 2)
#rangeθ = reshape([100.0], 1, 1)
rangeθ = reshape([220.0], 1, 1)
@info "rangeθ, rangel:" rangeθ, [lmin lmax]
# rangeθ = [0.111 25]
# build btg model

function load_unrate_btg()
    return btg(trainingData0, rangeθ, rangeλ; priorθ = inverseUniform(rangeθ), quadtype = myquadtype)
end

function getseed()
    return randseed
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module Matrix

using JuMP
using Cbc

include("Format.jl")
using .Format

include("FC.jl")
using .FC

include("Half.jl")
using .Half

export solve

# Matrix method solver -- uses CBC MIP optimizer
# Author: Jason Liu
function solve(m, s)
    model = Model(Cbc.Optimizer)
    set_optimizer_attribute(model, "logLevel", 0)
    set_optimizer_attribute(model, "allowableGap", 1e-5)
    set_optimizer_attribute(model, "cuts", "off")

    upper = min(fc(m, s, output = 0), half(m, s, output = 0))
    @variable(model, 0 <= x[i = 1:m, j = 1:s] <= 1)
    @variable(model, 0 <= y[i = 1:m, j = 1:s] <= 1, Int)
    @variable(model, 1 // 3 <= z <= upper, start = upper)

    @objective(model, Max, z)

    for i = 1:m, j = 1:s
        @constraint(model, x[i, j] + y[i, j] <= 1)
        @constraint(model, x[i, j] + y[i, j] >= 1 // s)
        @constraint(model, z <= x[i, j] + y[i, j])
    end

    for i = 1:m
        @constraint(model, sum(x[i, :]) == 1)
    end

    for j = 1:s
        @constraint(model, sum(x[:, j]) == m // s)
    end

    optimize!(model)

    return rationalize(value(z), tol = 1e-8)
end

end


text:
#GC related functions.
weaken_store(as::Vector) = begin
    for i = 1:length(as)
        as[i] = WeakRef(as[i])
    end
    nothing
end
strengthen_store(as::Vector) = begin
    filter!(x -> x != WeakRef(), as)
    for i = 1:length(as)
        as[i] = as[i].value
    end
    nothing
end
weaken_deps(sdgs::Vector) = begin
    for s in sdgs
        for i = 1:length(s.deps)
            s.deps[i] = WeakRef(s.deps[i])
        end
    end
    nothing
end
strengthen_deps(sdgs::Vector) = begin
    for s in sdgs
        filter!(x -> x != WeakRef(), s.deps)
        for i = 1:length(s.deps)
            if isa(s.deps[i], WeakRef)
                s.deps[i] = s.deps[i].value
            end
        end
    end
    nothing
end

weaken_getindexes() = begin
    for g in getindexes
        g.value = g.struct[map(value, g.args)...]
    end
    for g in getindexes
        if !isweak(g.struct)
            weaken(g.struct)
        end
    end
    nothing
end
strengthen_getindexes() = begin
    for g in getindexes
        if isweak(g.struct)
            strengthen(g.struct)
        end
    end
    for g in getindexes
        g.value = nothing
    end
    nothing
end

wr_value(wr::WeakRef) = wr.value
wr_value(sr) = sr

isdependent(parent::RV, child::Sample)   = parent == child.det
isdependent(parent::RV, child::Det)      = in(parent, child.args)
isdependent(parent::RV, child::GetIndex) =
    in(parent, child.args) || (parent == child.struct[map(value, child.args)...])

args(s::Sample) = {s.det}
args(d::Det) = {d.args...}
args(g::GetIndex) =
    {g.struct[map(value, g.args)...], g.args...}

condition_rec(parent, child::RV) = nothing
condition_rec(parent::RV, child::RV) = begin
    @assert isdependent(parent, child)
    i = findfirst(x -> x==wr_value(child), parent.deps)
    if isa(parent.deps[i], WeakRef)
        parent.deps[i] = parent.deps[i].value
        for g_parent in args(parent)
            condition_rec(g_parent, parent)
        end
    end
    nothing
end
condition_rec(cond::Sample) = begin
    condition_rec(cond.det, cond)
    nothing
end
condition_rec() = begin
    for cond in conditions
        condition_rec(cond)
    end
    nothing
end

apply_to_sdg(f::Function) = begin
    f(samples)
    f(dets)
    f(getindexes)
    f(conditions)
end
gc_church() = begin
    #Remove redundant dependencies.
    apply_to_sdg(remove_deps)
    gc_disable()
    #Weaken links in lazy structures.
    weaken_getindexes()
    #Weaken all dependent links.
    apply_to_sdg(weaken_deps)
    #Strengthen conditioned links.
    condition_rec()
    #Weaken links in samples, dets, etc.
    apply_to_sdg(weaken_store)
    gc_enable()
    gc()
    gc_disable()
    #Strengthen links in samples, dets, etc.
    apply_to_sdg(strengthen_store)
    apply_to_sdg(strengthen_deps)
    strengthen_getindexes()
    gc_enable()
end
n_samples() = length(samples)


text:
using DifferentialEquations
using Interpolations
using Plots

theta_res = 0.075;     # Residual water content
theta_sat = 0.287;     # Volumetric water content at saturation
vg_alpha = 0.027;      # Inverse of the air entry potential [/cm]
vg_n = 3.96;           # Pore-size distribution index
vg_m = 1;              # Exponent
Ksat = 34 / 3600;      # Hydraulic conductivity at saturation [cm/s]

params = [theta_res theta_sat vg_alpha vg_n vg_m Ksat];
# Van Genuchten Function

function vG(ψ; params=params)
# ----------------------------------
# van Genuchten [1980] relationships
# ----------------------------------
# --- Soil parameters
θ_res = params[1];   # Residual water content
θ_sat = params[2];   # Volumetric water content at saturation
α = params[3];       # Inverse of the air entry potential
n = params[4];           # Pore-size distribution index
m = params[5];           # Exponent
Ksat = params[6];        # Hydraulic conductivity at saturation

# --- Effective saturation [Se] for specified matric potential [psi]
if (ψ < 0)
   Se = (1 + (α * abs(ψ))^n)^-m
else()
   Se = 1
end

# --- Volumetric soil moisture [theta] for specified matric potential [psi]
θ = θ_res + (θ_sat - θ_res) * Se

# --- Hydraulic conductivity [K] for specified matric potential [psi]
if Se < 1
   #K = Ksat * sqrt(Se) * (1 - (1 - Se^(1/m))^m)^2
   K = Ksat * 1.175e6 / (1.175e6 + abs(ψ)^4.74);
else
   K = Ksat
end

# --- Specific moisture capacity [cap] for specified matric potential [psi]
if (ψ < 0)
   num = α * m * n * (θ_sat - θ_res) * (α * abs(ψ))^(n-1)
   den =  (1 + (α * abs(ψ))^n)^(m+1)
   cap = num / den
else
   cap = 0.0
end

return θ, K, cap
end

# Discretization (should be flexible in delta-z, need to try out later)
N = 150
# Let's do layers
z_bnd = zeros(N+1);
for i=2:length(z_bnd)
   z_bnd[i] = z_bnd[i-1]-i-1;
end
# Layer thickness:
Δz = z_bnd[1:end-1]-z_bnd[2:end];
# Layer center:
z = 0.5*(z_bnd[1:end-1]+z_bnd[2:end]);
# Half Layer plus/minus thickness
Δz_half = similar(z_bnd);
Δz_half[1]=z_bnd[1]-z[1];
for i = 2:N
   Δz_half[i] = z[i-1]-z[i];
end
Δz_half[N+1] = z[N]-z_bnd[N+1];

#Δx[1:10]=-1
# Set prior Ψ
Ψ_0 = zeros(N);
θ_0 = zeros(N);
θ_0[:] .= 0.1;

Ψ_0[:] .= -100;
#Ψ_0[50:end]=-
#Ψ_0[50:end] .= -

function matric_potential(θ::Number; params=params)
# --- Calculate psi for a given theta
   θ_res = params[1];    # Residual water content
   θ_sat = params[2];    # Volumetric water content at saturation
   α = params[3];        # Inverse of the air entry potential
   n = params[4];            # Pore-size distribution index
   m = params[5];            # Exponent
   if θ<=θ_res
      θ=max(θ_res,θ)+1e-10
   elseif θ>=θ_sat
      θ=min(θ_sat,θ)-1e-10
   end

   Se = (θ  - θ_res) / (θ_sat - θ_res);
   psi = -((Se^(-1/m) - 1)^(1/n)) / α;
end

# Still too many allocations here but we can fine-tune later
#This basically solves Eq 8.24 in Gordon's book
function soil_water(du,u,p,t)
   re = vG.(u, params=p)
   # need interpolation to get to half levels
   #theta = [x[1] for x in re]
   KK = [x[2] for x in re]
   C = [x[3] for x in re]
   #KK[end]=0.0
   for i=1:length(u)
      # Top layer
      if i==1
         Kp = 0.5*(KK[i]+KK[i+1])
         Km = 0.0037
         # Just as test, boundary condition of Psi=-20 here.
         du[i] = (Km/Δz_half[i]*(-20-u[i])-Kp/Δz_half[i+1]*(u[i]-u[i+1])+Km-Kp)/Δz[i]/C[i]
         # Just as test, set flux boundary here.

      # Bottom Layer
      elseif i==length(u)
         Km = 0.5*(KK[i]+KK[i-1])
         # Free drainage at the bottom (prefer to model groundwater later):
         #du[i] = (Km/Δz[i]*(u[i-1]-u[i])+Km-KK[i])/Δz[i]/C[i]
         du[i] = (Km/Δz[i]*(u[i-1]-u[i])+Km-KK[i])/Δz[i]/C[i]

      else
         Km = 0.5*(KK[i]+KK[i-1])
         Kp = 0.5*(KK[i]+KK[i+1])
         du[i] = (Km/Δz_half[i]*(u[i-1]-u[i])-Kp/Δz_half[i+1]*(u[i]-u[i+1])+Km-Kp)/Δz[i]/C[i]
      end
   end
end

tspan = (0.0,20*3600.0)
p = params
prob = ODEProblem(soil_water,Ψ_0,tspan,p)
alg = ImplicitEuler()
# Save every 10min
@time sol = solve(prob, alg,saveat=10*60, dtmax=5);

plot(sol[:,1:end], z/100)

max_stars_repo_name: ruillercoaster/LSM-SPAM

max_issues_repo_name: ruillercoaster/LSM-SPAM

max_forks_repo_name: ruillercoaster/LSM-SPAM


text: print("Hello World")

max_stars_repo_name: notVansh/Hello-World

max_issues_repo_name: notVansh/Hello-World

max_forks_repo_name: notVansh/Hello-World


text:
using CairoMakie, Colors, CSV, DataFrames

export polyelement, datadir, load_rocket_specs

const datadir = normpath(joinpath(@__DIR__, "../data"))

"""
    polyelement(cl::Color) -> PolyElement

This is useful to create legends.

    fig = Figure(resolution = (800, 600))
    elements = polyelement.(colors)
    Legend(fig[row, col], elements, labels, title)

Make sure that the elements match up with the labels.
This must be true `length(elements) == length(labels)`.
"""
function polyelement(color::Color)
    PolyElement(polycolor = color)
end

function load_rocket_specs()
    path = joinpath(datadir, "rocket-specifications-2021.csv")
    rocketspecs = CSV.File(path) |> DataFrame
end


text:
using Knockout, WebIO, Blink
using Test

using WebIO: node

cleanup = !AtomShell.isinstalled()
cleanup && AtomShell.install()

# write your own tests here
s = Observable(["a", "b", "c"])
t = node(:select, attributes = Dict("data-bind" => "options : options", "id" => "myselect"));
n = (knockout(t, ["options" => s]));

w = Window(Blink.@d(:show => false)); sleep(10.0)

body!(w, n); sleep(5.0)

@test Blink.@js(w, document.querySelector("#myselect").children.length) == 3
@test Blink.@js(w, document.querySelector("#myselect").children[0].value) == "a"
@test Blink.@js(w, document.querySelector("#myselect").children[1].value) == "b"
@test Blink.@js(w, document.querySelector("#myselect").children[2].value) == "c"

s[] = ["c", "d"]
sleep(1.0)

@test Blink.@js(w, document.querySelector("#myselect").children.length) == 2
@test Blink.@js(w, document.querySelector("#myselect").children[0].value) == "c"
@test Blink.@js(w, document.querySelector("#myselect").children[1].value) == "d"

cleanup && AtomShell.uninstall()


text:
# Test for valid tensor expressions and parse them.
const prime = Symbol("'")

# test for assignment (copy into existing tensor) or definition (create new tensor)
isassignment(ex::Expr) = ex.head == :(=) || ex.head == :(+=) || ex.head == :(-=)
function isdefinition(ex::Expr)
    #TODO: remove when := is removed
    # if ex.head == :(:=)
    #     warn(":= will likely be deprecated as assignment operator in Julia, use ≔ (\\coloneq + TAB) or go to http://github.com/Jutho/TensorOperations.jl to suggest ASCII alternatives", once=true, key=:warnaboutcoloneq)
    # end
    return ex.head == :(:=) || ex.head == :(≔)
end

# get left and right hand side from an expression or definition
function getlhsrhs(ex::Expr)
    if ex.head in (:(=), :(+=), :(-=), :(:=), :(≔))
        return ex.args[1], ex.args[2]
    else
        error("invalid assignment or definition $ex")
    end
end

# test for a valid index
function isindex(ex)
    if isa(ex, Symbol) || isa(ex, Int)
        return true
    elseif isa(ex, Expr) && ex.head == prime && length(ex.args) == 1
        return isindex(ex.args[1])
    else
        return false
    end
end

# test for a simple tensor object indexed by valid indices
istensor(ex) = false
function istensor(ex::Expr)
    if ex.head == :ref || ex.head == :typed_hcat
        return all(isindex, ex.args[2:end])
    elseif ex.head == :typed_vcat
        length(ex.args) == 3 || return false
        if isa(ex.args[2], Expr) && ex.args[2].head == :row
            all(isindex, ex.args[2].args) || return false
        else
            isindex(ex.args[2]) || return false
        end
        if isa(ex.args[3], Expr) && ex.args[3].head == :row
            all(isindex, ex.args[3].args) || return false
        else
            isindex(ex.args[3]) || return false
        end
        return true
    end
    return false
end

# test for a generalized tensor, i.e. with scalar multiplication and conjugation
isgeneraltensor(ex) = false
function isgeneraltensor(ex::Expr)
    if istensor(ex)
        return true
    elseif ex.head == :call && ex.args[1] == :+ && length(ex.args) == 2
        # unary plus
        return isgeneraltensor(ex.args[2])
    elseif ex.head == :call && ex.args[1] == :- && length(ex.args) == 2
        # unary minus
        return isgeneraltensor(ex.args[2])
    elseif ex.head == :call && ex.args[1] == :conj && length(ex.args) == 2
        # conjugation
        return isgeneraltensor(ex.args[2])
    elseif ex.head == :call && ex.args[1] == :adjoint && length(ex.args) == 2
        # adjoint
        return isgeneraltensor(ex.args[2])
    elseif ex.head == prime && length(ex.args) == 1
        # adjoint
        return isgeneraltensor(ex.args[1])
    elseif ex.head == :call && ex.args[1] == :transpose && length(ex.args) == 2
        # conjugation
        return isgeneraltensor(ex.args[2])
    elseif ex.head == :call && ex.args[1] == :* && length(ex.args) == 3
        # scalar multiplication
        if isscalarexpr(ex.args[2]) && isgeneraltensor(ex.args[3])
            return true
        elseif isscalarexpr(ex.args[3]) && isgeneraltensor(ex.args[2])
            return true
        end
    elseif ex.head == :call && ex.args[1] == :/ && length(ex.args) == 3
        # scalar multiplication
        if isscalarexpr(ex.args[3]) && isgeneraltensor(ex.args[2])
            return true
        end
    elseif ex.head == :call && ex.args[1] == :\ && length(ex.args) == 3
        # scalar multiplication
        if isscalarexpr(ex.args[2]) && isgeneraltensor(ex.args[3])
            return true
        end
    end
    return false
end

function hastraceindices(ex)
    obj,leftind, rightind, = makegeneraltensor(ex)
    allind = vcat(leftind, rightind)
    return length(allind) != length(unique(allind))
end

# test for a scalar expression, i.e. no indices
function isscalarexpr(ex::Expr)
    if ex.head == :call && ex.args[1] == :scalar
        return true
    elseif ex.head in (:ref, :typed_vcat, :typed_hcat)
        return false
    else
        return all(isscalarexpr, ex.args)
    end
end
isscalarexpr(ex::Symbol) = true
isscalarexpr(ex::Number) = true
isscalarexpr(ex) = false

# test for a tensor expression, i.e. something that can be evaluated to a tensor
function istensorexpr(ex)
    if isgeneraltensor(ex)
        return true
    elseif isa(ex, Expr) && ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :-)
        return all(istensorexpr, ex.args[2:end]) # all arguments should be tensor expressions (we are not checking matching indices yet)
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :*
        return any(istensorexpr, ex.args[2:end]) # at least one argument should be a tensor expression
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :/ && length(ex.args) == 3
        return istensorexpr(ex.args[2]) && isscalarexpr(ex.args[3])
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :\ && length(ex.args) == 3
        return istensorexpr(ex.args[3]) && isscalarexpr(ex.args[2])
    end
    return false
end

# convert an expression into a valid index
function makeindex(ex)
    if isa(ex, Symbol) || isa(ex, Int)
        return ex
    elseif isa(ex, Expr) && ex.head == prime && length(ex.args) == 1
        return Symbol(makeindex(ex.args[1]), "′")
    else
        error("not a valid index: $ex")
    end
end

# extract the tensor object itself, as well as its left and right indices
function maketensor(ex)
    if isa(ex, Expr) && (ex.head == :ref || ex.head == :typed_hcat)
        object = esc(ex.args[1])
        leftind = Any[makeindex(x) for x in ex.args[2:end]]
        rightind = Any[]
        return (object, leftind, rightind)
    elseif isa(ex, Expr) && ex.head == :typed_vcat
        length(ex.args) <= 3 || throw(ArgumentError("invalid tensor index expression: $ex"))
        object = esc(ex.args[1])
        if isa(ex.args[2], Expr) && ex.args[2].head == :row
            leftind = Any[makeindex(x) for x in ex.args[2].args]
        elseif ex.args[2] == :_
            leftind = Any[]
        else
            leftind = Any[makeindex(ex.args[2])]
        end
        if length(ex.args) > 2 && isa(ex.args[3], Expr) && ex.args[3].head == :row
            rightind = Any[makeindex(x) for x in ex.args[3].args]
        elseif length(ex.args) == 2 || ex.args[3] == :_
            rightind = Any[]
        else
            rightind = Any[makeindex(ex.args[3])]
        end
        return (object, leftind, rightind)
    end
    throw(ArgumentError())
end

# extract the tensor object itself, as well as its left and right indices, its scalar factor and its conjugation flag
function makegeneraltensor(ex)
    if istensor(ex)
        object, leftind, rightind = maketensor(ex)
        return (object, leftind, rightind, 1, false)
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :+ && length(ex.args) == 2 # unary plus: pass on
        return makegeneraltensor(ex.args[2])
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :- && length(ex.args) == 2 # unary minus: flip scalar factor
        (object, leftind, rightind, α, conj) = makegeneraltensor(ex.args[2])
        return (object, leftind, rightind, Expr(:call, :-, α), conj)
    elseif isa(ex, Expr) && ex.head == :call && ex.args[1] == :conj && length(ex.args) == 2 # conjugation: flip conjugation flag and conjugate scalar factor
        (object, leftind, rightind, α, conj) = makegeneraltensor(ex.args[2])
        return (object, leftind, rightind, Expr(:call, :conj, α), !conj)
    elseif ex.head == :call && ex.args[1] == :* && length(ex.args) == 3 # scalar multiplication: muliply scalar factors
        if isscalarexpr(ex.args[2]) && isgeneraltensor(ex.args[3])
            (object, leftind, rightind, α, conj) = makegeneraltensor(ex.args[3])
            return (object, leftind, rightind, Expr(:call, :*, makescalar(ex.args[2]), α), conj)
        elseif isscalarexpr(ex.args[3]) && isgeneraltensor(ex.args[2])
            (object, leftind, rightind, α, conj) = makegeneraltensor(ex.args[2])
            return (object, leftind, rightind, Expr(:call, :*, α, makescalar(ex.args[3])), conj)
        end
    elseif ex.head == :call && ex.args[1] == :/ && length(ex.args) == 3 # scalar multiplication: muliply scalar factors
        if isscalarexpr(ex.args[3]) && isgeneraltensor(ex.args[2])
            (object, leftind, rightind, α, conj) = makegeneraltensor(ex.args[2])
            return (object, leftind, rightind, Expr(:call, :/, α, makescalar(ex.args[3])), conj)
        end
    elseif ex.head == :call && ex.args[1] == :\ && length(ex.args) == 3 # scalar multiplication: muliply scalar factors
        if isscalarexpr(ex.args[2]) && isgeneraltensor(ex.args[3])
            (object, leftind, rightind, α, conj) = makegeneraltensor(ex.args[3])
            return (object, leftind, rightind, Expr(:call, :\, makescalar(ex.args[2]), α), conj)
        end
    end
    throw(ArgumentError("not a valid generalized tensor expression $ex"))
end

function makescalar(ex::Expr)
    if ex.head == :call && ex.args[1] == :scalar
        @assert length(ex.args) == 2 && istensorexpr(ex.args[2])
        return :(scalar($(deindexify(nothing, 0, ex.args[2], 1, [], [], true))))
    elseif ex.head == :call
        return Expr(ex.head, ex.args[1], map(makescalar, ex.args[2:end])...)
    else
        return Expr(ex.head, map(makescalar, ex.args)...)
    end
end
makescalar(ex::Symbol) = esc(ex)
makescalar(ex) = ex

# for any tensor expression, get the list of uncontracted indices that would remain after evaluating that expression
function getindices(ex::Expr)
    if istensor(ex)
        _,leftind,rightind = maketensor(ex)
        return unique2(vcat(leftind, rightind))
    elseif ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :-)
        return getindices(ex.args[2]) # getindices on any of the args[2:end] should yield the same result
    elseif ex.head == :call && ex.args[1] == :*
        indices = getindices(ex.args[2])
        for k = 3:length(ex.args)
            append!(indices, getindices(ex.args[k]))
        end
        return unique2(indices)
    elseif ex.head == :call && ex.args[1] == :/
        indices = getindices(ex.args[2])
    elseif ex.head == :call && ex.args[1] == :\
        indices = getindices(ex.args[3])
    elseif ex.head == :call && length(ex.args) == 2
        return getindices(ex.args[2])
    else
        return Vector{Any}()
    end
end
getindices(ex) = Vector{Any}()

# get all unique indices appearing in a tensor expression
function getallindices(ex::Expr)
    if istensor(ex)
        _,leftind,rightind = maketensor(ex)
        return unique(vcat(leftind, rightind))
    elseif !isempty(ex.args)
        return unique(mapreduce(getallindices, vcat, ex.args))
    else
        return Vector{Any}()
    end
end
getallindices(ex) = Vector{Any}()

function geteltype(ex::Expr)
    if istensor(ex)
        obj,_,_ = maketensor(ex)
        return Expr(:call, :eltype, obj)
    elseif ex.head == :call && (ex.args[1] == :+ || ex.args[1] == :- || ex.args[1] == :* || ex.args[1] == :/)
        if length(ex.args) > 2
            return Expr(:call, :promote_type, map(geteltype, ex.args[2:end])...)
        else
            return geteltype(ex.args[2])
        end
    elseif ex.head == :call && ex.args[1] == :conj
        return geteltype(ex.args[2])
    elseif isscalarexpr(ex)
        return :(typeof($ex))
    else
        # return :(eltype($ex)) # would probably lead to doing the same operation twice
        throw(ArgumentError("unable to determine eltype"))
    end
end
geteltype(ex) = Expr(:call,:typeof, esc(ex))


text:
import ..SamplingConstraints:
    SequentialSamplingConstraint,
    OrderedSamplingAlgorithm
import ..UncertainDatasets:
    AbstractUncertainValueDataset

"""
    resample(udata::AbstractUncertainValueDataset,
        sequential_constraint::SequentialSamplingConstraint;
        quantiles = [0.0001, 0.9999])

Resample a dataset by imposing a sequential sampling constraint.

Before drawing the realization, all furnishing distributions are truncated to the provided
`quantiles` range. This is to avoid problems in case some distributions have infinite
support.

"""
resample(udata::AbstractUncertainValueDataset,
    sequential_constraint::SequentialSamplingConstraint;
    quantiles = [0.0001, 0.9999])

"""
    resample(udata::AbstractUncertainValueDataset,
        constraint::Union{SamplingConstraint, Vector{SamplingConstraint}},
        sequential_constraint::SequentialSamplingConstraint;
        quantiles = [0.0001, 0.9999])

Resample a dataset by first imposing regular sampling constraints on the furnishing
distributions, then applying a sequential sampling constraint.

Before drawing the realization, all furnishing distributions are truncated to the provided
`quantiles` range. This is to avoid problems in case some distributions have infinite
support.
"""
resample(udata::AbstractUncertainValueDataset,
    constraint::Union{SamplingConstraint, Vector{SamplingConstraint}},
    sequential_constraint::SequentialSamplingConstraint;
    quantiles = [0.0001, 0.9999])


text:
module RocketConcatMapToOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: concat_map_to()" begin

    println("Testing: operator concat_map_to()")

    struct DummyType end

    @testset begin
        @test_throws InvalidSubscribableTraitUsageError of(0) |> concat_map_to(DummyType())
    end

    run_testset([
        (
            source      = from([ 0, 0 ]) |> concat_map_to(from([ 1, 2 ])),
            values      = @ts([ 1, 2, 1, 2, c ]),
            source_type = Int
        ),
        (
            source      = from(1:5) |> concat_map_to(of(0.0)),
            values      = @ts([ 0.0, 0.0, 0.0, 0.0, 0.0, c ]),
            source_type = Float64
        ),
        (
            source      = from(1:5) |> concat_map_to(faulted(Float64, "err")),
            values      = @ts(e("err")),
            source_type = Float64
        ),
        (
            source      = from(1:5) |> concat_map_to(completed(Float64)),
            values      = @ts(c),
            source_type = Float64
        ),
        (
            source      = from(1:5) |> concat_map_to(never(Float64)),
            values      = @ts(),
            source_type = Float64
        ),
        (
            source      = completed() |> concat_map_to(of(0)),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = faulted(String, "e") |> concat_map_to(of(0)),
            values      = @ts(e("e")),
            source_type = Int
        ),
        (
            source      = never() |> concat_map_to(of(0)),
            values      = @ts(),
            source_type = Int
        ),
        (
            source      = from([ 0, 0 ]) |> async(0) |> concat_map_to(from([ 1, 2 ])),
            values      = @ts([ 1, 2 ] ~ [ 1, 2 ] ~ c),
            source_type = Int
        ),
        (
            source      = from([ 0, 0 ]) |> async(0) |> concat_map_to(from([ 1, 2 ]) |> async(0)),
            values      = @ts([ 1 ] ~ [ 2 ] ~ [ 1 ] ~  [ 2 ] ~ c),
            source_type = Int
        ),
        (
            source      = from([ of(1), completed(Int), of(2) ]) |> concat_map_to(of(0)),
            values      = @ts([ 0, 0, 0, c ]),
            source_type = Int
        ),
        (
            source      = from([ of(1), completed(Int), of(2) ]) |> concat_map(Int) |> concat_map_to(of(0)),
            values      = @ts([ 0, 0, c ]),
            source_type = Int
        ),
        (
            source      = from([ of(1), faulted(Int, "err"), of(2) ]) |> concat_map_to(of(0)),
            values      = @ts([ 0, 0, 0, c ]),
            source_type = Int
        ),
        (
            source      = from([ of(1), faulted(Int, "err"), of(2) ]) |> concat_map(Int) |> concat_map_to(of(0)),
            values      = @ts([ 0, e("err") ]),
            source_type = Int
        )
    ])

end

end


text:
using SparseArrayKit
using Test, TestExtras, LinearAlgebra, Random
import TensorOperations

#=
generate a whole bunch of random contractions, compare with the dense result
=#
function randn_sparse(T::Type{<:Number}, sz::Dims, p = 0.5)
    a = SparseArray{T}(undef, sz)
    for I in keys(a)
        if rand() < p
            a[I] = randn(T)
        end
    end
    return a
end

@timedtestset "Basic linear algebra" begin
    MAX_DIM = 20;
    MAX_LEGS = 4;

    dims = ntuple(l->rand(1:MAX_DIM), MAX_LEGS)
    ar = randn_sparse(Float64, dims)
    ac = randn_sparse(ComplexF64, dims)

    slice = (1:rand(1:dims[1]), rand(1:dims[2]), rand(1:dims[3]):dims[3], rand(1:dims[4]))
    @test ar[slice...] == Array(ar)[slice...]

    α = randn(ComplexF64)
    β = randn(Float64)
    γ = 2
    @test @constinferred(α*ar) == α*Array(ar)
    @test @constinferred(β*ar) == β*Array(ar)
    @test @constinferred(γ*ar) == γ*Array(ar)
    @test @constinferred(ac*β) == Array(ac)*β
    @test @constinferred(ar + ac) == Array(ar) + Array(ac)
    @test @constinferred(ac - ar) == Array(ac) - Array(ar)
    mar = @constinferred(-ar)
    @test iszero(ar + mar)
    ar2 = copy(ar)
    ar2[first(nonzero_keys(ar2))] = 0
    @test !(ar2 == ar)
    @test @constinferred(zero(ar)) + ac == ac
    @test ac == convert(SparseArray, Array(ac))

    @test norm(ar + @constinferred(α*ar)) ≈ norm(ar)*abs(1+α)
    @test @constinferred(norm(ac + ac*α)) ≈ norm(ac)*abs(1+α)
    @test ar*α ≈ α*ar
    @test norm(ar*α) ≈ norm(ar)*abs(α);
    @test ac/β ≈ β\ac;
    @test norm(ac/γ) ≈ norm(ac)/abs(γ);
    @test norm(ac) ≈ sqrt(real(@constinferred(dot(ac, ac))))
    @test @constinferred(dot(ac, ar)) ≈ dot(Array(ac), Array(ar))
    @test dot(ar, ac) ≈ conj(dot(ac, ar))

    @test @constinferred(lmul!(β, copy(ar))) == β * ar == β * Array(ar)
    @test @constinferred(rmul!(copy(ac), α)) == ac * α == Array(ac) * α
    @test @constinferred(rmul!(copy(ac), β)) == ac * β == Array(ac) * β
    @test @constinferred(ldiv!(β, copy(ar))) == ldiv!(β, Array(ar)) ≈ β \ ar
    @test @constinferred(rdiv!(copy(ac), α)) == rdiv!(Array(ac), α) ≈ ac / α
    @test @constinferred(rdiv!(copy(ac), β)) == rdiv!(Array(ac), β) ≈ ac / β
    @test @constinferred(conj!(copy(ac))) == conj!(Array(ac))
    @test_throws InexactError rmul!(copy(ar), α)
    @test_throws InexactError ldiv!(α, copy(ar))

    p = randperm(MAX_LEGS)
    @test permutedims(ar, p) == permutedims(Array(ar), p)

    @test @constinferred(axpy!(α, copy(ar), copy(ac))) == α * ar + ac
    @test @constinferred(axpby!(α, copy(ar), β, copy(ac))) == α * ar + β * ac
end

@timedtestset "Basic matrix algebra" begin
    using SparseArrays
    N = 100
    a = sprandn(ComplexF64, N, N, 0.1)
    b = sprandn(N, N, 0.1)
    aa = @constinferred(SparseArray(a))
    bb = SparseArray(b)
    @test aa == a
    @test @constinferred(aa * bb) ≈ SparseArray(a * b)
    @test aa * bb ≈ Array(aa) * Array(bb)
    @test adjoint(aa) == adjoint(Array(aa))
    @test transpose(aa) == transpose(Array(aa))
    @test aa'*bb ≈ Array(aa)'*Array(bb)
    @test aa'*bb' ≈ Array(aa)'*Array(bb)'
    @test aa*bb' ≈ Array(aa)*Array(bb)'
    @test @constinferred(one(aa)) == one(Array(aa))
    @test norm(one(aa)) ≈ sqrt(N)
    @test one(aa) + zero(aa) == one(bb)
    @test adjoint!(copy(aa), bb) == SparseArray(adjoint(bb)) == adjoint(Array(bb))
    @test transpose!(copy(aa), bb) == SparseArray(transpose(bb)) == transpose(Array(bb))
end

@timedtestset "random contractions with eltype = $T" for T in (Float64, ComplexF64)

    MAX_CONTRACTED_INDICES = 10
    MAX_OPEN_INDICES = 5
    MAX_DIM = 5
    MAX_IND_PER_TENS = 3
    NUM_TESTS = 100

    for i in 1:NUM_TESTS
        contracted_indices = repeat(collect(1:rand(1:MAX_CONTRACTED_INDICES)),2)
        open_indices = collect(1:rand(1:MAX_OPEN_INDICES))
        dimensions = [repeat(rand(1:MAX_DIM,Int(length(contracted_indices)/2)),2);rand(1:MAX_DIM,length(open_indices))]

        #generate a random tensor network contraction
        tensors = SparseArray[]
        indices = Vector{Int64}[]
        conjlist = Bool[]
        while !isempty(contracted_indices) || !isempty(open_indices)
            num_inds = rand(1:min(MAX_IND_PER_TENS,length(contracted_indices)+length(open_indices)))

            cur_inds = Int64[]
            cur_dims = Int64[]

            for ind in 1:num_inds
                curind_index = rand(1:(length(contracted_indices)+length(open_indices)))

                if curind_index <= length(contracted_indices)
                    push!(cur_inds,contracted_indices[curind_index])
                    push!(cur_dims,dimensions[curind_index])
                    deleteat!(contracted_indices,curind_index)
                    deleteat!(dimensions,curind_index)
                else
                    tind = curind_index - length(contracted_indices)
                    push!(cur_inds,-open_indices[tind])
                    push!(cur_dims,dimensions[curind_index])
                    deleteat!(open_indices,tind)
                    deleteat!(dimensions,curind_index)
                end
            end

            push!(tensors, randn_sparse(T, tuple(cur_dims...)))
            push!(indices, cur_inds)
            push!(conjlist, rand([true,false]))
        end

        length(tensors) == 1 && continue # very rare - but possible

        sparse_result = TensorOperations.ncon(tensors,indices,conjlist)
        dense_result = TensorOperations.ncon(Array.(tensors),indices,conjlist)

        @test Array(sparse_result) ≈ dense_result
    end
end


text:
using Mustache
using DrWatson

jobname(args...) = replace(savename(args...;connector="_"),"="=>"_")
driverdir(args...) = normpath(projectdir("..",args...))

function jobdict(params)
  nr = params[:nr]
  nnodes = params[:nnodes]
  tpn = params[:tpn]
  dt = params[:dt]
  write_solution = params[:write_solution]
  write_solution_freq = params[:write_solution_freq]
  k = params[:k]
  degree = params[:degree]
  mrelax = params[:mrelax]
  nstep = params[:nstep]
  d=Dict(
  "q" => "normal",
  "o" => datadir(jobname(params,"o.txt")),
  "e" => datadir(jobname(params,"e.txt")),
  "walltime" => "04:00:00",
  "ncpus" => nnodes*48,
  "tpn" => tpn,
  "mem" => "$(nnodes*48*4)gb",
  "name" => jobname(params),
  "numrefs" => params[:numrefs],
  "n" => nnodes*tpn,
  "nthreads" => params[:nthreads],
  "nnodes" => nnodes,
  "dt" => dt,
  "ws" => write_solution,
  "write_solution_freq" => write_solution_freq,
  "k" => k,
  "nr" => nr,
  "degree" => degree,
  "mrelax" => mrelax,
  "nstep" => nstep,
  "projectdir" => driverdir(),
  "modules" => driverdir("modules.sh"),
  "title" => datadir(jobname(params)),
  "sysimage" => driverdir("GalewskyShallowWaterThetaMethod.so"),
  )
  println(d)
  d
end

allparams = Dict(
 :nnodes => [32,40,56,64,80,88],
 :numrefs => 10,
 :write_solution => false,
 :write_solution_freq => 4,
 :dt => 480.0,
 :k => 1,
 :degree => 4,
 :mrelax => 5000,
 :nstep => 10,
 :nr => 5,
 :nthreads => 1,
 :tpn=>48,
 )

template = read(projectdir("jobtemplate.sh"),String)

dicts = dict_list(allparams)

for params in dicts
  jobfile = datadir(jobname(params,"sh"))
  open(jobfile,"w") do io
    render(io,template,jobdict(params))
  end
end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
#=
Created on 20/01/2021 10:58:10
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Here, we illustrate type-based dispatch, using the type system
=#

abstract type ODESolver end

struct Euler <: ODESolver end
struct SecondOrder <: ODESolver end

function solve_ode(f′, (t₀, tₑ), y₀, Δt, method::ODESolver)
    tsteps = t₀:Δt:tₑ
    # initialize solution vector
    y = zeros(length(tsteps))
    # set initial value
    y[1] = y₀
    for (i, t) in enumerate(tsteps)
        i == 1 && continue
        yᵢ₋₁ = y[i-1]
        yᵢ = step(f′, t, yᵢ₋₁, Δt, method)
        y[i] = yᵢ
    end
    return y
end

step(f′, t, yᵢ₋₁, Δt, ::Euler) = yᵢ₋₁ + Δt * f′(t, yᵢ₋₁)

function step(f′, t, yᵢ₋₁, Δt, ::SecondOrder)
    y = yᵢ₋₁ + 0.5Δt * f′(t, yᵢ₋₁)  # half step
    return yᵢ₋₁ + Δt * f′(t, y)
end

f′(t, y) = -y + sin(t * y)

y_euler = solve_ode(f′, (0, 10), 1.0, 0.1, Euler())
y_second_order = solve_ode(f′, (0, 10), 1.0, 0.1, SecondOrder())

plot(y_euler, label="Euler", xlabel="t", ylabel="y")
plot!(y_second_order, label="second order")


text:
@testset "Discretization" begin
  @testset "FIST" begin
    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(2,1),(2,2),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == [5]

    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(2,1),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == [4]

    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == [2,4]

    𝒫 = Chain(P2[(0,0),(1,1),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == []

    points = P2[(0,0),(1,0),(1,1),(2,1),(2,2),(1,2),(0,0)]
    connec = connect.([(4,5,6),(3,4,6),(3,6,1),(1,2,3)], Triangle)
    target = UnstructuredMesh(points[1:end-1], connec)
    𝒫 = PolyArea(points)
    mesh = discretize(𝒫, FIST())
    @test mesh == target
  end
end


text:
# ---
# title: 1313. Decompress Run-Length Encoded List
# id: problem1313
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/decompress-run-length-encoded-list/description/>
# hidden: true
# ---
#
# We are given a list `nums` of integers representing a list compressed with
# run-length encoding.
#
# Consider each adjacent pair of elements `[freq, val] = [nums[2*i],
# nums[2*i+1]]` (with `i >= 0`).  For each such pair, there are `freq` elements
# with value `val` concatenated in a sublist. Concatenate all the sublists from
# left to right to generate the decompressed list.
#
# Return the decompressed list.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3,4]
#     Output: [2,4,4,4]
#     Explanation: The first pair [1,2] means we have freq = 1 and val = 2 so we generate the array [2].
#     The second pair [3,4] means we have freq = 3 and val = 4 so we generate [4,4,4].
#     At the end the concatenation [2] + [4,4,4] is [2,4,4,4].
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,1,2,3]
#     Output: [1,3,3]
#
#
#
#
# **Constraints:**
#
#   * `2 <= nums.length <= 100`
#   * `nums.length % 2 == 0`
#   * `1 <= nums[i] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# ==============================================================================
# Filtered Oscillator
#
# See: https://flowstar.org/benchmarks/filtered-oscillator/
# ==============================================================================

using HybridSystems, MathematicalSystems, MathematicalPredicates, LazySets,
      LinearAlgebra, Reachability, Polyhedra, Optim

import LazySets.HalfSpace
import LazySets.Approximations: overapproximate, OctDirections

"""
    function filtered_oscillator(n0::Int=4, time_horizon::Float64=99.,
                                 one_loop_iteration::Bool=false)

Construct the filtered-oscillator model.

### Input

- `n0`                 -- (optional; default: `4`) number of filters
- `time_horizon`       -- (optional; default: `99.`) time horizon
- `one_loop_iteration` -- (optional; default: `false`) option to enforce a
                          single loop (see Notes section below)

### Notes

The option `one_loop_iteration` can be used to add an additional variable `k`.
This variable starts at value `1` and is multiplied by `2` when taking the
transition from location `ℓ₃` to location `ℓ₄`.
In location `ℓ₄` we add the invariant constraint `k ≤ 2.1`.
Thus location `ℓ₄` can only be entered once.
"""
function filtered_oscillator(n0::Int=4,
                             time_horizon::Float64=99.,
                             one_loop_iteration::Bool=false)

    n1 = (one_loop_iteration ? n0 + 1 : n0)
    n = n1 + 2
    z = zeros(n1)

    # transition graph (automaton)
    a = LightAutomaton(4)
    add_transition!(a, 3, 4, 1)
    add_transition!(a, 4, 2, 2)
    add_transition!(a, 2, 1, 3)
    add_transition!(a, 1, 3, 4)

    # common flow
    A = zeros(n, n)
    A[1,1], A[2,2] = -2., -1.
    A[3,1], A[3,3] = 5., -5.
    for i = 4 : n-1
        A[i,i-1], A[i,i] = 5., -5.
    end

    # modes

    # mode 1
    b = [1.4; -0.7; z]
    X = HPolyhedron([HalfSpace([-0.714286; -1.0; z], 0.0),  # 0.714286*x + y >= 0
                     HalfSpace([1.0; 0.0; z], 0.0)])  # x <= 0
    m_1 = CACS(A, b, X)

    # mode 2
    b = [-1.4; 0.7; z]
    X = HPolyhedron([HalfSpace([1.0; 0.0; z], 0.0),  # x <= 0
                     HalfSpace([0.714286; 1.0; z], 0.0)])  # 0.714286*x + y <= 0
    m_2 = CACS(A, b, X)

    # mode 3
    b = [1.4; -0.7; z]
    X = HPolyhedron([HalfSpace([-1.0; 0.0; z], 0.0),  # x >= 0
                     HalfSpace([-0.714286; -1.0; z], 0.0)])  # 0.714286*x + y >= 0
    m_3 = CACS(A, b, X)

    # mode 4
    b = [-1.4; 0.7; z]
    X = HPolyhedron([HalfSpace([0.714286; 1.0; z], 0.0),  # 0.714286*x + y <= 0
                     HalfSpace([-1.0; 0.0; z], 0.0)])  # x >= 0
    if one_loop_iteration
        # k <= 2 (2.1 to avoid numerical issues)
        addconstraint!(X, HalfSpace([zeros(n-1); 1.], 2.1))
    end
    m_4 = CACS(A, b, X)

    m = [m_1, m_2, m_3, m_4]

    # transitions

    # transition l3 -> l4
    X_l3l4 = HPolyhedron([HalfSpace([-1.0; 0.0; z], 0.0),  # x >= 0
                          HalfSpace([-0.714286; -1.0; z], 0.0),  # 0.714286*x + y >= 0
                          HalfSpace([0.714286; 1.0; z], 0.0)])  # 0.714286*x + y <= 0
    if one_loop_iteration
        A_trans_34 = Matrix(1.0I, n, n)
        A_trans_34[n, n] = 2.  # k' = k * 2
        r1 = ConstrainedLinearMap(A_trans_34, X_l3l4)
    else
        r1 = ConstrainedIdentityMap(n, X_l3l4)
    end

    # transition l4 -> l2
    X_l4l2 = HPolyhedron([HalfSpace([0.714286; 1.0; z], 0.0),  # 0.714286*x + y <= 0
                          HalfSpace([-1.0; 0.0; z], 0.0),  # x >= 0
                          HalfSpace([1.0; 0.0; z], 0.0)])  # x <= 0
    r2 = ConstrainedIdentityMap(n, X_l4l2)

    # transition l2 -> l1
    X_l2l1 = HPolyhedron([HalfSpace([1.0; 0.0; z], 0.0),  # x <= 0
                          HalfSpace([-0.714286; -1.0; z], 0.0),  # 0.714286*x + y >= 0
                          HalfSpace([0.714286; 1.0; z], 0.0)])  # 0.714286*x + y <= 0
    r3 = ConstrainedIdentityMap(n, X_l2l1)

    # transition l1 -> l3
    X_l1l3 = HPolyhedron([HalfSpace([-0.714286; -1.0; z], 0.0),  # 0.714286*x + y >= 0
                          HalfSpace([-1.0; 0.0; z], 0.0),  # x >= 0
                          HalfSpace([1.0; 0.0; z], 0.0)])  # x <= 0
    r4 = ConstrainedIdentityMap(n, X_l1l3)

    r = [r1, r2, r3, r4]

    # switchings
    s = [HybridSystems.AutonomousSwitching()]

    HS = HybridSystem(a, m, r, s)

    # initial condition in mode 3
    low = [0.2; -0.1; zeros(n1)]
    high = [0.3; 0.1; zeros(n1)]
    if one_loop_iteration
        low[end] = 1.0
        high[end] = 1.0
    end
    X0 = Hyperrectangle(low=low, high=high)

    problem = InitialValueProblem(HS, [(3, X0)])

    # safety property
    border = HPolyhedron([HalfSpace([0.0; -1.0; z], -0.5)])  # y >= 0.5
    property = is_disjoint_from(border)

    options = Options(:T=>time_horizon, :mode=>"reach", :property=>property,
                      :verbosity=>0)

    solver_options = Options(:vars=>1:n, :δ=>0.01, :plot_vars=>[1, 2],
                             :ε_proj=>0.001, :project_reachset=>false)

    return (problem, options, solver_options)
end

"""
    get_projection(sol, projected_dims)

Overapproximate and project a flowpipe.

### Input

- `sol`            -- a flowpipe solution
- `projected_dims` -- an integer vector of length 2 with the variables to project

### Output

A `ReachSolution` containing the overapproximated and projected set.

### Notes

The overapproximation used is octagonal directions.
"""
function get_projection(sol, projected_dims)
    n = sol.options[:n] # system's dimension
    oa = x -> overapproximate(x, OctDirections(n))
    sol_oa = ReachSolution([ReachSet(CartesianProductArray([oa(rs.X)]), rs.t_start, rs.t_end) for rs in sol.Xk], sol.options)
    sol_proj = ReachSolution(project_reach(sol_oa.Xk, projected_dims, n, sol.options), sol.options)
    return sol_proj
end

# single run of filtered oscillator
problem, options, solver_options = filtered_oscillator()
result = solve(problem, options, BFFPSV18(solver_options), ApproximatingDiscretePost());

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
module SimulatedNeuralMoments
# the type that holds the model specifics
struct SNMmodel
    modelname::String # name of model
    lb::Vector # vector of lower bounds. Can be -Inf, if desired
    ub::Vector # vector of upper bounds. Can be inf, if desired
    insupport::Function # function that checks if a draw is valid
    prior::Function # function that evaluates the prior at draw
    priordraw::Function # function that returns a draw from prior
    auxstat::Function # function that returns an array of draws of statistic
end

include("MakeNeuralMoments.jl")
include("FromEconometrics.jl")
include("SNM.jl")
include("MCMC.jl")
export SNMmodel, MakeNeuralMoments
export prettyprint, dstats
export TransformStats, NeuralMoments, EstimateΣ, H
export MCMC
end


text:
module TrajTools

export Trajectory
export dist_tmp, matmul_tmp

include("read_trajec.jl")
include("distances.jl")
include("calcs.jl")

const dist_tmp = zeros(MVector{3})
const matmul_tmp = zeros(MVector{3})

end


text:
using Test
using ManifoldsBase
using ManifoldsBase: AbstractNumbers, ℝ, ℂ

struct DummyPowerRepresentation <: AbstractPowerRepresentation end

@testset "Power Manifold" begin
    M = ManifoldsBase.DefaultManifold(3)
    N = PowerManifold(M, NestedPowerRepresentation(), 2)
    p = [zeros(3), ones(3)]
    q = [ones(3), zeros(3)]
    @testset "Constructors" begin
        @test repr(N) ==
              "PowerManifold(DefaultManifold(3; field = ℝ), NestedPowerRepresentation(), 2)"
        # add to type
        @test repr(N^3) ==
              "PowerManifold(DefaultManifold(3; field = ℝ), NestedPowerRepresentation(), 2, 3)"
        # add to type
        @test repr(PowerManifold(N, 3)) ==
              "PowerManifold(DefaultManifold(3; field = ℝ), NestedPowerRepresentation(), 2, 3)"
        # switch type
        @test repr(PowerManifold(N, DummyPowerRepresentation(), 3)) ==
              "PowerManifold(DefaultManifold(3; field = ℝ), DummyPowerRepresentation(), 2, 3)"
        # nest
        @test repr(PowerManifold(N, NestedPowerRepresentation(), 3)) ==
              "PowerManifold(PowerManifold(DefaultManifold(3; field = ℝ), NestedPowerRepresentation(), 2), NestedPowerRepresentation(), 3)"
    end
    @test N^3 == PowerManifold(M, NestedPowerRepresentation(), 2, 3)
    @test ManifoldsBase.get_iterator(N^3) == Base.product(Base.OneTo(2), Base.OneTo(3))

    @testset "point/tangent checks" begin
        pE1 = [zeros(3), ones(2)] # one component wrong
        pE2 = [zeros(2), ones(2)] # both wrong
        @test is_manifold_point(N, p)
        @test is_manifold_point(N, p, true)
        @test !is_manifold_point(N, pE1)
        @test !is_manifold_point(N, pE2)
        @test_throws ComponentManifoldError is_manifold_point(N, pE1, true)
        @test_throws CompositeManifoldError is_manifold_point(N, pE2, true)
        # tangent - test base
        @test is_tangent_vector(N, p, p)
        @test is_tangent_vector(N, p, p, true)
        @test !is_tangent_vector(N, pE1, p)
        @test !is_tangent_vector(N, pE2, p)
        # tangents - with proper base
        @test is_tangent_vector(N, p, p, true)
        @test !is_tangent_vector(N, p, pE1)
        @test !is_tangent_vector(N, p, pE2)
        @test_throws ComponentManifoldError is_tangent_vector(N, p, pE1, true)
        @test_throws CompositeManifoldError is_tangent_vector(N, p, pE2, true)
    end
    @testset "specific functions" begin
        @test distance(N, p, q) == sqrt(sum(distance.(Ref(M), p, q) .^ 2))
        @test exp(N, p, q) == p .+ q
        @test retract(N, p, q) == p .+ q
        @test ManifoldsBase.get_iterator(N) == Base.OneTo(2)
        @test injectivity_radius(N) == injectivity_radius(M)
        @test injectivity_radius(N, p) == injectivity_radius(M, p)
        p2 = allocate(p)
        copyto!(N, p2, p)
        @test !(p === p2)
        @test p == p2
        q2 = allocate(q)
        copyto!(N, q2, p, q)
        @test !(q === q2)
        @test q == q2
        @test inner(N, p, q, q) == sum(inner.(Ref(M), p, q, q))
        @test isapprox(N, p, q) == (all(isapprox.(Ref(M), p, q)))
        @test isapprox(N, p, p) == (all(isapprox.(Ref(M), p, p)))
        @test isapprox(N, p, q, p) == (all(isapprox.(Ref(M), p, q, p)))
        @test isapprox(N, p, p, p) == (all(isapprox.(Ref(M), p, p, p)))
        @test log(N, p, q) == q .- p
        @test inverse_retract(N, p, q) == q .- p
        @test manifold_dimension(N) == 2 * manifold_dimension(M)
        @test mid_point(N, p, q) == mid_point.(Ref(M), p, q)
        @test sqrt(inner(N, p, q, q)) == norm(N, p, q)
        @test project(N, p) == p
        @test project(N, p, q) == q
        @test power_dimensions(N) == (2,)
        @test power_dimensions(N^3) == (2, 3)
        m = ParallelTransport()
        @test vector_transport_to(N, p, p, q) == p
        @test vector_transport_to(N, p, p, q, m) == p
        @test vector_transport_to(N, p, p, q, PowerVectorTransport(m)) == p
        q2 = [zeros(3), zeros(3)]
        vector_transport_to!(N, q2, p, q, p)
        @test q2 == q
        @test vector_transport_direction(N, p, p, q) == p
        @test vector_transport_direction(N, p, p, q, m) == p
        @test vector_transport_direction(N, p, p, q, PowerVectorTransport(m)) == p
        q2 = [zeros(3), zeros(3)]
        vector_transport_direction!(N, q2, p, q, p)
        @test q2 == q
        @test p[N, 1] == p[1]
        p[N, 1] = 2 .* ones(3)
        @test p[N, 1] == 2 .* ones(3)
        @test p[N, [2, 1]] == [p[2], p[1]]
        @test view(p, N, 2) isa SubArray
        @test view(p, N, 2) == p[2]
    end
    @testset "Basis, coordinates & vector" begin
        v = get_coordinates(N, p, q, DefaultBasis())
        @test v == [q[1]..., q[2]...]
        @test get_vector(N, p, v, DefaultBasis()) == q
        B = get_basis(N, p, DefaultBasis())
        @test get_coordinates(N, p, q, B) == v
        v2 = zeros(size(v))
        get_coordinates!(N, v2, p, q, B)
        @test v2 == v
        @test get_coordinates(N, p, q, DefaultOrthonormalBasis()) == v
        @test B.data.bases[1].data == get_basis(M, p[1], DefaultBasis()).data
        @test B.data.bases[2].data == get_basis(M, p[2], DefaultBasis()).data
        @test get_vector(N, p, v, B) == q
        @test zero_tangent_vector(N, p) == [zeros(3), zeros(3)]
        B2 = DiagonalizingOrthonormalBasis([ones(3), ones(3)])
        B3 = get_basis(N, p, B2)
        @test sprint(show, "text/plain", B) == """$(DefaultBasis()) for a power manifold
        Basis for component (1,):
        $(sprint(show, "text/plain", B.data.bases[1]))
        Basis for component (2,):
        $(sprint(show, "text/plain", B.data.bases[2]))
        """
    end
    @testset "Zero index manifold" begin
        M = ManifoldsBase.DefaultManifold()
        N = PowerManifold(M, NestedPowerRepresentation(), 3)
        p = [1.0, 2.0, 3.0]

        @test p[N, 1] == 1.0
        @test zero_tangent_vector(N, p) == zero(p)
    end
end


text:
function save_flux_model(filepath::String, p::AbstractVector, s::AbstractVector, re)
    base_name, ext = splitext(filepath)
    if ext == ".jld2"
        save(filepath, Dict("model" => re(p, s) |> cpu))
        @info "Saved model to " * filepath
        return nothing
    end
    dict = Dict("parameters" => p |> cpu, "states" => s |> cpu)
    ps_filepath = base_name * ".ps" * ext
    serialize(ps_filepath, dict)
    @info "Model Parameters and States saved to " * ps_filepath
    re_filepath = base_name * ".re" * ext
    serialize(re_filepath, re)
    @info "Restructure Function saved to " * re_filepath
    return nothing
end

save_flux_model(filepath::String, m) = save_flux_model(filepath, destructure_parameters_states(m)...)

function load_flux_model(filepath::String, re=nothing)
    base_name, ext = splitext(filepath)
    if ext == ".jld2"
        model = load(filepath, "model")
        @info "Model loaded from " * filepath
        re === nothing && return model
        p, s, _ = destructure_parameters_states(model)
        return re(p, s)
    end
    ps_filepath = base_name * ".ps" * ext
    re_filepath = base_name * ".re" * ext
    ps = deserialize(ps_filepath)
    parameters, states = ps["parameters"], ps["states"]
    if re === nothing
        try
            re = deserialize(re_filepath)
        catch e
            @warn "Unable to load the restructure function from " * re_filepath * ". Please pass it manually."
            rethrow(e)
        end
    end
    return re(parameters, states)
end

# mutable struct ModelCheckpointer
#     checkpoint_frequency::Int
#     step::Int
#     checkpoint_filepath::String
#     maximum_checkpoints::Int
#     checkpoint_counter::Int

# end