text:
# simulation step
function tick!(
    state::Tuple{AbstractGraph, AbstractArray}, tweet_list::AbstractArray,
    tick_nr::Int64, config::Config
    )

    agent_list = state[2]
    for agent_idx in shuffle(1:length(agent_list))
        this_agent = agent_list[agent_idx]
        if this_agent.active && (rand() < this_agent.check_regularity)
            update_feed!(state, agent_idx, config)
            update_perceiv_publ_opinion!(state, agent_idx)
            update_opinion!(state, agent_idx, config)
            like(state, agent_idx, config)
            retweet!(state, agent_idx, config)
            drop_input!(state, agent_idx, config)
            add_input!(state, agent_idx, config)
            inclin_interact = deepcopy(this_agent.inclin_interact)
            while inclin_interact > 0
                if rand() < inclin_interact
                    publish_tweet!(state, tweet_list, agent_idx, tick_nr)
                end
                inclin_interact -= 1.0
            end
            update_check_regularity!(state, agent_idx, config)
            this_agent.inactive_ticks = 0
        elseif this_agent.active
            this_agent.inactive_ticks += 1
            if this_agent.inactive_ticks > config.simulation.max_inactive_ticks
                set_inactive!(state, agent_idx, tweet_list)
            end
        end
    end
    update_network!(state, config)
    return log_network(state, tick_nr)
end

# the actual simulation
function simulate(config::Config = Config())
    graph = create_network(config.network.agent_count,config.network.m0)
    init_state = (graph, create_agents(graph))
    state = deepcopy(init_state)
    tweet_list = Array{Tweet, 1}(undef, 0)
    graph_list = Array{AbstractGraph, 1}([graph])
    df = DataFrame(
        TickNr = Int64[],
        AgentID = Int64[],
        Opinion = Float64[],
        PerceivPublOpinion = Float64[],
        InclinInteract = Float64[],
        InactiveTicks = Int64[],
        Indegree = Float64[],
        ActiveState = Bool[]
    )
    for i in 1:config.simulation.n_iter
        append!(df, tick!(state, tweet_list, i, config))
        if i % ceil(config.simulation.n_iter / 10) == 0
            print(".")
            push!(graph_list, deepcopy(state[1]))
        end
    end

    tweet_df = DataFrame(
                Opinion = [t.opinion for t in tweet_list],
                Weight = [t.weight for t in tweet_list],
                Source_Agent = [t.source_agent for t in tweet_list],
                Published_At = [t.published_at for t in tweet_list],
                Likes = [t.like_count for t in tweet_list],
                Retweets = [t.retweet_count for t in tweet_list]
                )

    return (df, tweet_df, graph_list), state, init_state
end
;

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
using StatsBase

lags = 1:40
# ! A bit of overhead here...
AC = FeatureSet([x->autocor(x, [ℓ]; demean=true)[1]::Float64 for ℓ ∈ lags],
                Symbol.(["AC_$ℓ" for ℓ ∈ lags]),
                [["correlation"] for ℓ ∈ lags],
                ["Autocorrelation at lag $ℓ" for ℓ ∈ lags])
export AC

# ! Lots of overhead here; the pacf function has to recalculate for all smaller lags.
# ! Time to think about SuperFeatures...
AC_partial = FeatureSet([x->pacf(x, [ℓ]; method=:regression)[1]::Float64 for ℓ ∈ lags],
                Symbol.(["AC_partial_$ℓ" for ℓ ∈ lags]),
                [["correlation"] for ℓ ∈ lags],
                ["Partial autocorrelation at lag $ℓ (regression method)" for ℓ ∈ lags])
export AC_partial

ext: jl


text:
typealias AbstractTensorSample AbstractSample{SecondOrder}
#############################################################################################################################################################
### TensorSample is used by samplers that compute (up to) the tensor of the log-target (for ex SmMALA)

immutable TensorSample{N<:Number,T<:AbstractFloat,V<:AbstractVector,A<:AbstractArray,D<:AbstractArray} <: AbstractTensorSample
    values::A
    loglikelihood::V
    logprior::V
    gradloglikelihood::A
    gradlogprior::A
    tensorloglikelihood::D
    tensorlogprior::D
    TensorSample(v::AbstractArray{N},ll::AbstractVector{T},lp::AbstractVector{T},gl::AbstractArray{T},gp::AbstractArray{T},
                 tl::AbstractArray{T},tp::AbstractArray{T}) = new(v,ll,lp,gl,gp,tl,tp)
end

function TensorSample{N<:Number,T<:AbstractFloat}(nparas::Integer,nsamples::Integer,::Type{N},::Type{T})
    v = zeros(N,_valuestuple(nparas,nsamples))
    u = _tensortuple(nparas,nparas,nsamples)
    TensorSample{N,T,Vector,Array,Array}(v,Vector{T}(nsamples),Vector{T}(nsamples),similar(v,T),similar(v,T),similar(v,T,u),similar(v,T,u))
end

### Factory function
_samples{N<:Number,T<:AbstractFloat}(::Type{Val{:tensor}},nparas::Integer,nsamples::Integer,::Type{N},::Type{T}) = TensorSample(nparas,nsamples,N,T)

@inline sampletypename(s::TensorSample) = :tensor

################################################################################################################################################################
### TangentTensorSample is used by samplers that approximate the calculation of the tensor of the log-target (for ex TrSmMALARandom)

immutable TangentTensorSample{N<:Number,T<:AbstractFloat,V<:AbstractVector,A<:AbstractArray,D<:AbstractArray} <: AbstractTensorSample
    values::A
    loglikelihood::V
    logprior::V
    gradloglikelihood::A
    gradlogprior::A
    tensorloglikelihood::D
    tensorlogprior::D
    tangentvectors::D
    TangentTensorSample(v::AbstractArray{N},ll::AbstractVector{T},lp::AbstractVector{T},gl::AbstractArray{T},gp::AbstractArray{T},
                        tl::AbstractArray{T},tp::AbstractArray{T},tv::AbstractArray{T}) = new(v,ll,lp,gl,gp,tl,tp,tv)
end

function TangentTensorSample{N<:Number,T<:AbstractFloat}(nparas::Integer,nsamples::Integer,ntangents::Integer,::Type{N},::Type{T})
    v = zeros(N,_valuestuple(nparas,nsamples))
    u1 = _tensortuple(nparas,nparas,nsamples)
    u2 = _tensortuple(nparas,ntangents,nsamples)
    TangentTensorSample{N,T,Vector,Array,Array}(v,Vector{T}(nsamples),Vector{T}(nsamples),similar(v,T),similar(v,T),similar(v,T,u1),similar(v,T,u1),similar(v,T,u2))
end

_samples{N<:Number,T<:AbstractFloat}(::Type{Val{:tangent}},nparas::Integer,nsamples::Integer,::Type{N},::Type{T},ntangents::Integer) = TangentTensorSample(nparas,nsamples,ntangents,N,T)

### Additional functionality for tangent tensor samples
numtangents(s::TangentTensorSample) = size(s.tangentvectors,2)

similar(s::TangentTensorSample,nsamples::Integer =numsamples(s),tosamplename::Symbol =sampletypename(s)) = (args=tosamplename==:tangent?(numtangents(s),):() ; _samples(Val{tosamplename},numparas(s),nsamples,sampletype(s),calculationtype(s),args...))

@inline sampletypename(s::TangentTensorSample) = :tangent


text:
using Documenter, OutputCollectors

makedocs(
    modules = [OutputCollectors],
    sitename = "OutputCollectors.jl",
)

deploydocs(
    repo = "github.com/JuliaPackaging/OutputCollectors.jl.git",
    push_preview = true,
)


text:
## Julia test script to run on Travis
## workaround until Pkg3 can REQUIRE master versions of packages

## for these packages, use master
Pkg.clone("https://github.com/tpapp/ContinuousTransformations.jl.git");

if VERSION ≥ v"0.6-"
    ## workaround for https://github.com/JuliaPlots/Plots.jl/issues/753
    for pkg in ["Plots", "PlotThemes", "PlotUtils"]
        Pkg.add(pkg)
        Pkg.checkout(pkg)
    end
end

## original test script
Pkg.clone(pwd()); Pkg.build("ParametricFunctions");
Pkg.test("ParametricFunctions"; coverage=true)


text:
using HTTP, Gumbo

const PAGE_URL = "https://en.wikipedia.org/wiki/Julia_(programming_language)"
const LINKS = String[]

function fetchpage(url)
  response = HTTP.get(url)
  if response.status == 200 && parse(Int, Dict(response.headers)["Content-Length"]) > 0
    String(response.body)
  else
    ""
  end
end

function extractlinks(elem)
  if  isa(elem, HTMLElement) &&
      tag(elem) == :a &&
      in("href", collect(keys(attrs(elem))))
        url = getattr(elem, "href")
        startswith(url, "/wiki/") && ! occursin(":", url) && push!(LINKS, url)
  end

  for child in children(elem)
    extractlinks(child)
  end
end

content = fetchpage(PAGE_URL)

if ! isempty(content)
  dom = Gumbo.parsehtml(content)

  extractlinks(dom.root)
end

display(unique(LINKS))

max_stars_repo_name: PacktPublishing/Julia-Projects

max_issues_repo_name: PacktPublishing/Julia-Projects

max_forks_repo_name: PacktPublishing/Julia-Projects


text:
# ---
# title: 689. Maximum Sum of 3 Non-Overlapping Subarrays
# id: problem689
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, Dynamic Programming
# link: <https://leetcode.com/problems/maximum-sum-of-3-non-overlapping-subarrays/description/>
# hidden: true
# ---
#
# In a given array `nums` of positive integers, find three non-overlapping
# subarrays with maximum sum.
#
# Each subarray will be of size `k`, and we want to maximize the sum of all
# `3*k` entries.
#
# Return the result as a list of indices representing the starting position of
# each interval (0-indexed). If there are multiple answers, return the
# lexicographically smallest one.
#
# **Example:**
#
#
#
#     Input: [1,2,1,2,6,7,5,1], 2
#     Output: [0, 3, 5]
#     Explanation: Subarrays [1, 2], [2, 6], [7, 5] correspond to the starting indices [0, 3, 5].
#     We could have also taken [2, 1], but an answer of [1, 3, 5] would be lexicographically larger.
#
#
#
#
# **Note:**
#
#   * `nums.length` will be between 1 and 20000.
#   * `nums[i]` will be between 1 and 65535.
#   * `k` will be between 1 and floor(nums.length / 3).
#
#
#
#
## @lc code=start
using LeetCode

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


text:
## Tests involing building packages and whatnot

@testset "Builder Packaging" begin
    # Gotta set this guy up beforehand
    tarball_path = nothing
    tarball_hash = nothing

	# Test building with both `make` and `cmake`
    for script in (libfoo_make_script, libfoo_cmake_script)
        product_storage = tempname()
        mkpath(product_storage)

		begin
			build_path = tempname()
			mkpath(build_path)
			prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], platform)
			cd(joinpath(dirname(@__FILE__),"build_tests","libfoo")) do
				run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)
				@test build(ur, "foo", libfoo_products(prefix), script, platform, prefix)
			end

			# Next, package it up as a .tar.gz file
            tarball_path, tarball_hash = package(prefix, joinpath(product_storage, "libfoo"), v"1.0.0"; verbose=true)
			@test isfile(tarball_path)

			# Delete the build path
			rm(build_path, recursive = true)
		end

		# Test that we can inspect the contents of the tarball
		contents = list_tarball_files(tarball_path)
		@test "bin/fooifier" in contents
		@test "lib/libfoo.$(Libdl.dlext)" in contents

		# Install it within a new Prefix
		temp_prefix() do prefix
			# Install the thing
			@test install(tarball_path, tarball_hash; prefix=prefix, verbose=true)

			# Ensure we can use it
			fooifier_path = joinpath(bindir(prefix), "fooifier")
			libfoo_path = joinpath(libdir(prefix), "libfoo.$(Libdl.dlext)")
			check_foo(fooifier_path, libfoo_path)
		end

        rm(product_storage; recursive=true, force=true)
    end
end

if lowercase(get(ENV, "BINARYBUILDER_FULL_SHARD_TEST", "false") ) == "true"
    @info("Beginning full shard test...")
    # Perform a sanity test on each and every shard.
    @testset "Shard sanity tests" begin
        for shard_platform in expand_gcc_versions(supported_platforms())
            @info(" --> $(triplet(shard_platform))")
            build_path = tempname()
            mkpath(build_path)

            # build with make
            prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], shard_platform)
            cd(joinpath(dirname(@__FILE__),"build_tests","libfoo")) do
                run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

                # Build libfoo, warn if we fail
                @test build(ur, "foo", libfoo_products(prefix), libfoo_make_script, shard_platform, prefix)
            end
            rm(build_path, recursive = true)

            # build again with cmake
            mkpath(build_path)
            prefix, ur = BinaryBuilder.setup_workspace(build_path, [], [], [], shard_platform)
			cd(joinpath(dirname(@__FILE__),"build_tests","libfoo")) do
                run(`cp $(readdir()) $(joinpath(prefix.path,"..","srcdir"))/`)

                # Build libfoo, warn if we fail
                @test build(ur, "foo", libfoo_products(prefix), libfoo_cmake_script, shard_platform, prefix)
            end
            rm(build_path, recursive = true)
        end
    end
end

# Testset to make sure we can build_tarballs() from a local directory
@testset "build_tarballs() local directory based" begin
    build_path = tempname()
    local_dir_path = joinpath(build_path, "libfoo")
    mkpath(local_dir_path)

    cd(build_path) do
        # Just like we package up libfoo into a tarball above, we'll just copy it
        # into a new directory and use build_tarball's ability to auto-package
        # local directories to do all the heavy lifting.
        libfoo_dir = joinpath(@__DIR__, "build_tests", "libfoo")
        run(`cp -r $(libfoo_dir)/$(readdir(libfoo_dir)) $local_dir_path`)

        build_tarballs(
            [], # fake ARGS
            "libfoo",
            v"1.0.0",
            [local_dir_path],
            libfoo_make_script,
            [Linux(:x86_64, :glibc)],
            libfoo_products,
            [], # no dependencies
        )

        # Make sure that worked
        @test isfile("products/libfoo.v1.0.0.x86_64-linux-gnu.tar.gz")
        @test isfile("products/build_libfoo.v1.0.0.jl")
    end
end

# Testset to make sure we can build_tarballs() from a git repository
@testset "build_tarballs() Git-Based" begin
    # Skip this testset on Travis, because its libgit2 is broken right now.
    if get(ENV, "TRAVIS", "") == "true"
        return
    end

    build_path = tempname()
    git_path = joinpath(build_path, "libfoo.git")
    mkpath(git_path)

    cd(build_path) do
        # Just like we package up libfoo into a tarball above, we'll create a fake
        # git repo for it here, then build from that.
        repo = LibGit2.init(git_path)
        LibGit2.commit(repo, "Initial empty commit")
        libfoo_dir = joinpath(@__DIR__, "build_tests", "libfoo")
        run(`cp -r $(libfoo_dir)/$(readdir(libfoo_dir)) $git_path/`)
        for file in ["fooifier.cpp", "libfoo.c", "Makefile"]
            LibGit2.add!(repo, file)
        end
        commit = LibGit2.commit(repo, "Add libfoo files")

        # Now build that git repository for Linux x86_64
        sources = [
            git_path =>
            LibGit2.string(LibGit2.GitHash(commit)),
        ]

        build_tarballs(
            [], # fake ARGS
            "libfoo",
            v"1.0.0",
            sources,
            "cd libfoo\n$libfoo_make_script",
            [Linux(:x86_64, :glibc)],
            libfoo_products,
            [], # no dependencies
        )

        # Make sure that worked
        @test isfile("products/libfoo.v1.0.0.x86_64-linux-gnu.tar.gz")
        @test isfile("products/build_libfoo.v1.0.0.jl")
    end

    rm(build_path; force=true, recursive=true)
end

@testset "build_tarballs() --only-buildjl" begin
    build_path = tempname()
    mkpath(build_path)
    cd(build_path) do
        # Clone down OggBuilder.jl
        repo = LibGit2.clone("https://github.com/staticfloat/OggBuilder", ".")

        # Check out a known-good tag
        LibGit2.checkout!(repo, string(LibGit2.GitHash(LibGit2.GitCommit(repo, "v1.3.3-6"))))

        # Reconstruct binaries!  We don't want it to pick up BinaryBuilder.jl information from CI,
        # so wipe out those environment variables through withenv:
        blacklist = ["CI_REPO_OWNER", "CI_REPO_NAME", "TRAVIS_REPO_SLUG", "TRAVIS_TAG", "CI_COMMIT_TAG"]
        withenv((envvar => nothing for envvar in blacklist)...) do
            m = Module(:__anon__)
            Core.eval(m, quote
                ARGS = ["--only-buildjl"]
            end)
            Base.include(m, joinpath(build_path, "build_tarballs.jl"))
        end

        # Read in `products/build.jl` to get download_info
        m = Module(:__anon__)
        download_info = Core.eval(m, quote
            using BinaryProvider
            # Override BinaryProvider functionality so that it doesn't actually install anything
            function install(args...; kwargs...); end
            function write_deps_file(args...; kwargs...); end
        end)
        # Include build.jl file to extract download_info
        Base.include(m, joinpath(build_path, "products", "build_Ogg.v1.3.3.jl"))
        download_info = Core.eval(m, :(download_info))

        # Test that we get the info right about some of these platforms
        bin_prefix = "https://github.com/staticfloat/OggBuilder/releases/download/v1.3.3-6"
        @test download_info[Linux(:x86_64)] == (
            "$bin_prefix/Ogg.v1.3.3.x86_64-linux-gnu.tar.gz",
            "6ef771242553b96262d57b978358887a056034a3c630835c76062dca8b139ea6",
        )
        @test download_info[Windows(:i686)] == (
            "$bin_prefix/Ogg.v1.3.3.i686-w64-mingw32.tar.gz",
            "3f6f6f524137a178e9df7cb5ea5427de6694c2a44ef78f1491d22bd9c6c8a0e8",
        )
    end
end


text:
using Revise
using AdFem
using PyPlot
using LinearAlgebra
using ADCME
using Statistics
using SparseArrays

function get_fault_info(mmesh, c1, c2)
    # c1=(0.5, 0.0)
    # c2=(1.5, 0.5)
    vec = [c2[1]-c1[1], c2[2]-c1[2]]
    k = vec[2]/(vec[1]+1e-8)
    normal_vec = [-vec[2] vec[1]]./norm(vec)
    nodes = mmesh.nodes
    ids = Int64[]
    coords = Any[]
    dist(x0, y0) = abs(k * (x0-c1[1]) - (y0-c1[2])) / sqrt(k^2 + 1)
    for i = 1:size(nodes, 1)
        x, y = nodes[i,:]
        if dist(x, y) < 1e-3 && (c1[1] <= x <= c2[1])
            push!(ids, i)
            push!(coords, [x, y])
        end
    end
   return ids, hcat(coords...), normal_vec
end

"""
isredblack(elem_idx, node_idx, mmesh)::Bool
"""
function fix_mesh(mmesh, red_fault_idx, isredblack)
    nodes = mmesh.nodes
    elems = mmesh.elems

    extra_nodes = [nodes;nodes[red_fault_idx,:]]
    black_fault_idx =collect( mmesh.nnode .+ (1:length(red_fault_idx)) )
    black_fault_idx_map = Dict()
    for i = 1:length(red_fault_idx)
        black_fault_idx_map[red_fault_idx[i]] = black_fault_idx[i]
    end

    conn = zeros(Int64, length(red_fault_idx), 2)
    for i = 1:length(red_fault_idx)
        conn[i,:] = [red_fault_idx[i] black_fault_idx[i]]
    end

    is_fault_node = Set(red_fault_idx)
    for i = 1:mmesh.nelem
        candidate = elems[i,:]
        for j = 1:3
            if elems[i,j] in is_fault_node
                if isredblack(i, elems[i,j], mmesh)
                    continue
                else
                    candidate[j] = black_fault_idx_map[elems[i,j]]
                end
            end
        end
        elems[i,:] = candidate
    end
    Mesh(extra_nodes, elems), conn
end

function rectangle_redblack(elem_idx, node_idx, mmesh)
    node = mean(mmesh.nodes[mmesh.elems[elem_idx,:],:], dims=1)
    x = node[1,1]
    if x<1.0
        return true
    else
        return false
    end
end

function fault_redblack(elem_idx, node_idx, mmesh)
    node = mean(mmesh.nodes[mmesh.elems[elem_idx,:],:], dims=1)
    x, y = node
    return y-0.5x+0.25>0
end

function viscoelasticity_solver(invη, μ, λ, Forces, mmesh)

    ## alpha-scheme
    β = 1/4; γ = 1/2
    a = b = 0.1    # damping

    fn_G = invη->begin
      G = tensor([1/Δt+2/3*μ*invη -μ/3*invη 0.0
        -μ/3*invη 1/Δt+2/3*μ*invη 0.0
        0.0 0.0 1/Δt+μ*invη])
      invG = inv(G)
    end
    invG = map(fn_G, invη)
    S = tensor([2μ/Δt+λ/Δt λ/Δt 0.0
        λ/Δt 2μ/Δt+λ/Δt 0.0
        0.0 0.0 μ/Δt])
    H = invG*S

    M = compute_fem_mass_matrix1(mmesh)
    Zero = spzero(mmesh.ndof, mmesh.ndof)
    M = [M Zero;Zero M]

    K = compute_fem_stiffness_matrix(H, mmesh)
    C = a*M + b*K # damping matrix
    L = M + γ*Δt*C + β*Δt^2*K
    bddof = [bdnode; bdnode .+ mmesh.ndof]
    L, _ = impose_Dirichlet_boundary_conditions(L, zeros(2mmesh.ndof), bddof, zeros(length(bddof)))
    L = factorize(L)

    a = TensorArray(NT+1); a = write(a, 1, zeros(2mmesh.ndof))
    v = TensorArray(NT+1); v = write(v, 1, zeros(2mmesh.ndof))
    d = TensorArray(NT+1); d = write(d, 1, zeros(2mmesh.ndof))
    U = TensorArray(NT+1); U = write(U, 1, zeros(2mmesh.ndof))
    Sigma = TensorArray(NT+1); Sigma = write(Sigma, 1, zeros(get_ngauss(mmesh), 3))
    Varepsilon = TensorArray(NT+1); Varepsilon = write(Varepsilon, 1, zeros(get_ngauss(mmesh), 3))

    function condition(i, tas...)
      i <= NT
    end

    function body(i, tas...)
      a_, v_, d_, U_, Sigma_, Varepsilon_ = tas
      a = read(a_, i)
      v = read(v_, i)
      d = read(d_, i)
      U = read(U_, i)
      Sigma = read(Sigma_, i)
      Varepsilon = read(Varepsilon_, i)

      res = batch_matmul(invG/Δt, Sigma)
      F = compute_strain_energy_term(res, mmesh) - K * U
      rhs = Forces[i] - F

      td = d + Δt*v + Δt^2/2*(1-2β)*a
      tv = v + (1-γ)*Δt*a
      rhs = rhs - C*tv - K*td
      rhs = scatter_update(rhs, bddof, zeros(length(bddof)))

      ## alpha-scheme
      a = L\rhs # bottleneck
      d = td + β*Δt^2*a
      v = tv + γ*Δt*a
      U_new = d

      Varepsilon_new = eval_strain_on_gauss_pts(U_new, mmesh)
      Sigma_new = update_stress_viscosity(Varepsilon_new, Varepsilon, Sigma, invη, μ*ones(get_ngauss(mmesh)), λ*ones(get_ngauss(mmesh)), Δt)

      i+1, write(a_, i+1, a), write(v_, i+1, v), write(d_, i+1, d), write(U_, i+1, U_new),
            write(Sigma_, i+1, Sigma_new), write(Varepsilon_, i+1, Varepsilon_new)
    end

    i = constant(1, dtype=Int32)
    _, _, _, _, u, sigma, varepsilon = while_loop(condition, body,
                      [i, a, v, d, U, Sigma, Varepsilon])

    U = stack(u)
    Sigma = stack(sigma)
    Varepsilon = stack(varepsilon)

    U, Sigma, Varepsilon
end

function compute_slip_boundary_condition( conn, mmesh)
    ii = Int64[]
    jj = Int64[]
    vv = Float64[]
    k = 0
    # u - u' = du
    # v - v' = dv
    for i = 1:size(conn, 1)
        k += 1
        push!(ii, k)
        push!(jj, conn[i,1])
        push!(vv, 1.0)
        push!(ii, k)
        push!(jj, conn[i,2])
        push!(vv, -1.0)
    end

    for i = 1:size(conn, 1)
        k += 1
        push!(ii, k)
        push!(jj, conn[i,1] + mmesh.ndof)
        push!(vv, 1.0)
        push!(ii, k)
        push!(jj, conn[i,2] + mmesh.ndof)
        push!(vv, -1.0)
    end
    sparse(ii, jj, vv, 2size(conn,1), 2mmesh.ndof)
end

function subduction_solver(invη, gravity, Conn, slipvec, bdnode, mmesh)
    slipvec = constant(slipvec)
    Conn = constant(Conn)
    invη = constant(invη)
    @assert size(slipvec,1)==NT
    @assert size(slipvec,2)==size(Conn,1)
    @assert size(Conn,2) == 2mmesh.ndof
    @assert isa(gravity, Real)
    @assert size(invη)==(get_ngauss(mmesh),)

    l_length = size(Conn, 1)
    ## alpha-scheme
    β = 1/4; γ = 1/2
    a = b = 0.1    # damping

    fn_G = invη->begin
      G = tensor([1/Δt+2/3*μ*invη -μ/3*invη 0.0
        -μ/3*invη 1/Δt+2/3*μ*invη 0.0
        0.0 0.0 1/Δt+μ*invη])
      invG = inv(G)
    end
    invG = map(fn_G, invη)
    S = tensor([2μ/Δt+λ/Δt λ/Δt 0.0
        λ/Δt 2μ/Δt+λ/Δt 0.0
        0.0 0.0 μ/Δt])
    H = invG*S

    M = compute_fem_mass_matrix1(mmesh)
    Zero = spzero(mmesh.ndof, mmesh.ndof)
    M = [M Zero;Zero M]

    K = compute_fem_stiffness_matrix(H, mmesh)
    C = a*M + b*K # damping matrix
    L = M + γ*Δt*C + β*Δt^2*K
    L = [L Conn'
        Conn spzero(l_length, l_length)]
    bddof = [bdnode; bdnode .+ mmesh.ndof]
    L, _ = impose_Dirichlet_boundary_conditions(L, zeros(2mmesh.ndof+l_length), bddof, zeros(length(bddof)))
    L = factorize(L)

    a = TensorArray(NT+1); a = write(a, 1, zeros(2mmesh.ndof))
    v = TensorArray(NT+1); v = write(v, 1, zeros(2mmesh.ndof))
    d = TensorArray(NT+1); d = write(d, 1, zeros(2mmesh.ndof))
    U = TensorArray(NT+1); U = write(U, 1, zeros(2mmesh.ndof))
    Sigma = TensorArray(NT+1); Sigma = write(Sigma, 1, zeros(get_ngauss(mmesh), 3))
    Varepsilon = TensorArray(NT+1); Varepsilon = write(Varepsilon, 1, zeros(get_ngauss(mmesh), 3))

    Forces = compute_fem_source_term(zeros(get_ngauss(mmesh)), -ones(get_ngauss(mmesh))*gravity, mmesh)

    function condition(i, tas...)
      i <= NT
    end

    function body(i, tas...)
      a_, v_, d_, U_, Sigma_, Varepsilon_ = tas
      a = read(a_, i)
      v = read(v_, i)
      d = read(d_, i)
      U = read(U_, i)
      Sigma = read(Sigma_, i)
      Varepsilon = read(Varepsilon_, i)

      res = batch_matmul(invG/Δt, Sigma)
      F = compute_strain_energy_term(res, mmesh) - K * U
      rhs = Forces - F

      td = d + Δt*v + Δt^2/2*(1-2β)*a
      tv = v + (1-γ)*Δt*a
      rhs = rhs - C*tv - K*td
      rhs = scatter_update(rhs, bddof, zeros(length(bddof)))

      ## alpha-scheme
      rhs = [rhs; slipvec[i]]
      a = L\rhs # bottleneck
      a = a[1:2mmesh.ndof]
      d = td + β*Δt^2*a
      v = tv + γ*Δt*a
      U_new = d

      Varepsilon_new = eval_strain_on_gauss_pts(U_new, mmesh)
      Sigma_new = update_stress_viscosity(Varepsilon_new, Varepsilon, Sigma, invη, μ*ones(get_ngauss(mmesh)), λ*ones(get_ngauss(mmesh)), Δt)

      i+1, write(a_, i+1, a), write(v_, i+1, v), write(d_, i+1, d), write(U_, i+1, U_new),
            write(Sigma_, i+1, Sigma_new), write(Varepsilon_, i+1, Varepsilon_new)
    end

    i = constant(1, dtype=Int32)
    _, _, _, _, u, sigma, varepsilon = while_loop(condition, body,
                      [i, a, v, d, U, Sigma, Varepsilon])

    U = stack(u)
    Sigma = stack(sigma)
    Varepsilon = stack(varepsilon)

    U, Sigma, Varepsilon
end


text:
using CUDA
@inline function CheckBlockBounds(offset, block, grids)
        x::Int32 = offset.x + block.x
        if x < 1 || x > grids
                return true
        end
        x = offset.y + block.y
        if x < 1 || x > grids
                return true
        end
        x = offset.z + block.z
        if x < 1 || x > grids
                return true
        end
        return false
end
function cuGridKnnSimpleNoInbounds(Points, Queries,
        Pointsperblock, Queriesperblock,
        IntPointsperblock, IntQueriesperblock,
        Distances, Neighbours,
        dimensions, grids, offset)

    if CheckBlockBounds(offset, blockIdx(), grids)
        return nothing
    end
    q_bid::Int32 = ((blockIdx().z - 1)*gridDim().y*gridDim().x
            + (blockIdx().y - 1)*gridDim().x
            + blockIdx().x)
    p_bid::Int32 = ((blockIdx().z + offset.z- 1)*gridDim().y*gridDim().x
            + (blockIdx().y + offset.y - 1)*gridDim().x
            + blockIdx().x + offset.x)
    tid::Int32 = threadIdx().x + (threadIdx().y-1)*blockDim().x
    stride::Int32 = (blockDim().x)*(blockDim().y)

    SharedPoints = @cuDynamicSharedMem(Float32, (dimensions, stride))
    SharedQueries = @cuDynamicSharedMem(Float32, (dimensions, stride), dimensions*stride*sizeof(Float32))
    # query = @view SharedQueries[:, tid]

    # startPoints::Int32 = 0
    # (threadIdx().x == 1)  && ( startPoints = IntPointsperblock[p_bid])
    # startPoints = shfl_sync(0xffffffff, startPoints, 1);
     startPoints::Int32 = IntPointsperblock[p_bid]
    startQueries::Int32 = IntQueriesperblock[q_bid]
    totalPoints::Int32 = Pointsperblock[p_bid]
    totalQueries::Int32 = Queriesperblock[q_bid]
    dist::Float32 = 0
    nb::Int32 = 0
    # 32 regs
    for q::Int32 = 0:stride:(totalQueries - 1)
        ## Read queries
        qIndex::Int32 = startQueries + q + tid
        if tid + q <= totalQueries
            for d::Int32 = 1:dimensions
                    SharedQueries[d,tid] = Queries[qIndex, d]
            end
            dist = Distances[qIndex]
            nb = Neighbours[qIndex]
        end
        for p::Int32 = 0:stride:(totalPoints - 1)
            ## Read points
            sync_threads()
            if tid + p <= totalPoints
                for d::Int32 = 1:dimensions
                    SharedPoints[d, tid] = Points[ startPoints + p + tid, d]
                end
            end
            sync_threads()
            ## Calculate and compare distance
            bounds::Int32 = CUDA.min(stride, totalPoints-p)
            for i::Int32 = 1:bounds
                tempdist::Float32 = 0
                for d::Int32 = 1:dimensions
                    tempdist += CUDA.pow(SharedQueries[d,tid]-SharedPoints[d, (i+tid-2)%bounds+1],2)
                end
                tempdist = CUDA.sqrt(tempdist)
                if tempdist < dist
                    dist = tempdist
                    nb = startPoints + p + (i+tid-2)%bounds+1
                end
            end
        end
        if tid + q <= totalQueries
            Distances[qIndex] = dist
            Neighbours[qIndex] = nb
        end
    end
    return nothing
end

function cuda_knn_simple_no_inbounds(OrderedPoints, OrderedQueries,
    PointsPerBlock, QueriesPerBlock,
    IntegralPointsPerBlock, IntegralQueriesPerBlock,
    numOfPoints, numOfQueries, numOfGrids, dimensions)

    ## Data transfer
    devPoints = CuArray(OrderedPoints)
    devQueries = CuArray(OrderedQueries)
    devPointsPerBlock = CuArray(PointsPerBlock)
    devQueriesPerBlock = CuArray(QueriesPerBlock)
    devIntegralPointsPerBlock = CuArray(IntegralPointsPerBlock)
    devIntegralQueriesPerBlock = CuArray(IntegralQueriesPerBlock)
    devRes = CUDA.fill(Float32(100), numOfQueries)
    devNeighbours = CuArray(zeros(Int32, numOfQueries))
    ## Config
    dimensions = Int32(dimensions)
    grids = Int32(numOfGrids)
    offset = (x=Int32(0),y=Int32(0),z=Int32(0))
    thread_groups = 2
    shmem=2*thread_groups*32*dimensions*sizeof(Float32)
    ## Running for main block.
    @cuda(blocks=(grids,grids,grids), threads=(32,thread_groups), shmem=shmem,
      cuGridKnnSimpleNoInbounds(devPoints, devQueries,
        devPointsPerBlock, devQueriesPerBlock,
        devIntegralPointsPerBlock, devIntegralQueriesPerBlock,
        devRes, devNeighbours, dimensions, grids, offset))
    #Running for neighbour blocks.
    # for x::Int32 = -1:1
    #     for y::Int32 = -1:1
    #         for z::Int32 = -1:1
    #             if x == 0 && y == 0 && z == 0
    #                 continue
    #             end
    #             offset = (x=x,y=y,z=z)
    #             @cuda(blocks=(grids,grids,grids), threads=(32,thread_groups), shmem=shmem,
    #               cuGridKnnSimpleNoInbounds(devPoints, devQueries,
    #                 devPointsPerBlock, devQueriesPerBlock,
    #                 devIntegralPointsPerBlock, devIntegralQueriesPerBlock,
    #                 devRes, devNeighbours, dimensions, grids, offset))
    #         end
    #     end
    # end

    return Array(devNeighbours), Array(devRes)

end


text:
using Distributions, Gadfly, Compose, Cairo, DataFrames, StatsBase, Distances, Colors

function gen_colors(n) # to create your own colors, here based on one of the brewer series
    cs = distinguishable_colors(n,
        [colorant"#66c2a5", colorant"#fc8d62", colorant"#8da0cb", colorant"#e78ac3",
            colorant"#a6d854", colorant"#ffd92f", colorant"#e5c494", colorant"#b3b3b3"],
        lchoices=Float64[58, 45, 72.5, 90],
        transform=c->deuteranopic(c, 0.1),
        cchoices=Float64[20,40],
        hchoices=[75,51,35,120,180,210,270,310]
    )
    convert(Vector{Color}, cs)
end

n_agents = 50
n_steps = 50
start = rand(n_agents)

bc_ar = Array{Float64,2}(undef, n_agents, n_steps + 1)

function bc_upd(ϵ::Float64, α::Float64, τ::Float64, averaging)
    bc_ar[:, 1] = start
    for j in 1:n_steps, i in 1:n_agents
        bc_ar[i, j + 1] = α*averaging(bc_ar[:, j][findall(abs.(bc_ar[:, j] .- bc_ar[:, j][i]) .< ϵ)]) + (1 - α)*τ
    end
    return bc_ar
end

ϵ = 0.25
α = 0.9
τ = 0.7

res = bc_upd(ϵ, α, τ, mean)

df = DataFrame(res')
rename!(df, [Symbol("$i") for i in 1:n_agents])
df = stack(df)
df[!, :steps] = repeat(1:n_steps + 1, outer=n_agents)

p = plot(df, x=:steps, y=:value, color=:variable, Geom.point, Geom.line,
    Coord.cartesian(xmax=40),
    Guide.xlabel("Time"),
    Guide.ylabel("Opinion"),
    Guide.title("ϵ = $ϵ / α = $α / τ = $τ / arithmetic mean"),
    yintercept=[τ], Geom.hline(style=:dot, color=colorant"grey"),
    Guide.annotation(compose(context(), Compose.text(22, τ + 0.015, "τ"), fontsize(13pt))),
    Theme(key_position=:none, point_size=1.25pt))

draw(PNG("BCU0250907.png", 7inch, 5inch), p)


text:
"""
This is a sample `DataConsumer` that is used for ingesting the parsed data from `LazyCSV`.
It works in two modes:
  - untyped: in this mode, each record is read and its fields are parsed, but no type
             type checking is done on them and the string value of fields are printed to
             the output stream. The only error detected in this mode is having an incorrect
             number of fields in a record.
  - typed: in this mode, fields are parsed to their desired format and type errors are
           detected. Finally, the values are printed to the output stream in a table-like
           format, where the type of each column is also printed in the header.
"""
struct PrintConsumer{REC_TYPE <:RecordType, IO_TP <: IO} <: DataConsumer
    out::IO_TP
    rec_type::REC_TYPE
    headers::Vector{String}
end

const UntypedPrintConsumer{IO_TP} = PrintConsumer{UntypedRecord, IO_TP}
function UntypedPrintConsumer(out::IO_TP, num_fields::Int=-1) where {IO_TP}
    UntypedPrintConsumer{IO_TP}(out, UntypedRecord(num_fields), Vector{String}())
end

const TypedPrintConsumer{IO_TP} = PrintConsumer{TypedRecord, IO_TP}
function TypedPrintConsumer(out::IO_TP, field_types::Type{T},
                            fallback_type::FieldType=NO_TYPE,
                            other_type_handlers=Tuple{}()) where {IO_TP, T}
    TypedPrintConsumer{IO_TP}(out,
                              TypedRecord(field_types, fallback_type, other_type_handlers),
                              Vector{String}())
end

function PrintConsumer(out::IO_TP, num_fields::Int=-1) where {IO_TP}
    UntypedPrintConsumer(out, num_fields)
end

function consume_header(pc::PrintConsumer, f::File, iter_res)
    consume_header(pc, f, iter_res) do field_header, i
        push!(pc.headers, string(field_header))
    end
end

function produce_header(pc::UntypedPrintConsumer, f::File)
    i = 0
    for field_header in pc.headers
        i += 1
        csv_field_string(pc.out, f, field_header, i, true)
    end
    j = i
    for n in (i+1):length(record_length(pc.rec_type))
        j += 1
        csv_field_string(pc.out, f, "$n", n, true)
    end
    j > 0 && write(pc.out, "\n")
    pc
end

function type_str(tp::FieldType)
    types = Vector{String}()
    if tp & INT_TYPE      != 0; push!(types, "int"); end
    if tp & FLOAT_TYPE    != 0; push!(types, "float"); end
    if tp & DATE_TYPE     != 0; push!(types, "date"); end
    if tp & DATETIME_TYPE != 0; push!(types, "datetime"); end
    if tp & BOOL_TYPE     != 0; push!(types, "bool"); end
    if tp & CHAR_TYPE     != 0; push!(types, "char"); end
    if tp & STRING_TYPE   != 0; push!(types, "string"); end
    if tp & MISSING_TYPE  != 0; push!(types, "missing"); end
    join(types, ",")
end

type_str_in_header(tp::FieldType) = " ($(type_str(tp)))"

function type_str_in_header(pc::TypedPrintConsumer, i::Int)
    type_str_in_header(field_type(pc.rec_type, i))
end

function produce_header(pc::TypedPrintConsumer, f::File)
    i = 0
    for field_header in pc.headers
        i += 1
        csv_field_string(pc.out, f, "$(strip(field_header))$(type_str_in_header(pc, i))",
                         i, true)
    end

    for n in (i+1):length(record_length(pc.rec_type))
        csv_field_string(pc.out, f, "$n$(type_str_in_header(pc, n))", n, true)
    end
    write(pc.out, "\n")
    pc
end
function consume_field(pc::UntypedPrintConsumer, f::File, pos::FilePos,
                       line, field_str, index::Int)
    csv_field_string(pc.out, f, strip(field_str), index)
    index == length(first_rec_fields(f)) && write(pc.out, "\n")
    true
end

function consume_field(pc::TypedPrintConsumer, f::File, pos::FilePos,
                       line, field_str, index::Int)
    field_tp = field_type(pc.rec_type, index)
    force_quote = field_tp in [[STRING_TYPE], [CHAR_TYPE]]
    field_value = strip(field_str)
    (isvalid, valid_type) = has_valid_type(field_tp, field_value)
    if isvalid
        csv_field_string(pc.out, f, field_value, index, force_quote)
        index == length(first_rec_fields(f)) && write(pc.out, "\n")
    end
    isvalid
end

function consume_field_error(pc::PrintConsumer, f::File, pos::FilePos,
                             line, field_str, index::Int)
    csv_field_string(pc.out, f, "FIELD_ERR($(strip(field_str)))", index, true)
    index == length(first_rec_fields(f)) && write(pc.out, "\n")
end

function consume_rec_error(pc::PrintConsumer, f::File, pos::FilePos, line_str)
    write(pc.out, "ERROR >> ")
    if SHOW_DETAILED_ERROR
        write(pc.out, "$(length(first_rec_fields(f))) fields:")
        for (i, field) in enumerate(first_rec_fields(f))
            write(pc.out, "\n    $i ->")
            write(pc.out, field)
        end
        write(pc.out, " in \nline: ")
    end
    write(pc.out, string(line_str))
end

function csv_string(buff::IO, csv_file::File, consumer::DataConsumer = PrintConsumer(buff))
    consume(consumer, csv_file)
    for pos_line in csv_file
        (pos, line) = pos_line
        i = 1
        for field in first_rec_fields(csv_file)
            csv_field_string(buff, csv_file, field, i)
            i += 1
        end
        write(buff, "\n")
    end
end

"""
Parses a CSV file in an untyped way and prints it back to the output stream.
If a row has an inconsistent number of fields, it will be reported as an error.
"""
function csv_string(csv_file::File, num_fields::Int=-1)
    buff = IOBuffer()
    csv_string(buff, csv_file, PrintConsumer(buff, num_fields))
    seekstart(buff)
    read(buff, String)
end

"""
Parses a CSV file in an typed way (by parsing the values to their intended types) and prints
it back to the output stream.
If a row has an inconsistent number of fields, it will be reported as an error.
Also, if a field has an icorrect type, it'll also be reported as an error.
"""
function typed_csv_string(csv_file::File, field_types::Type{T},
                          fallback_type::FieldType=NO_TYPE,
                          other_type_handlers=Tuple{}()) where {T}
    buff = IOBuffer()
    csv_string(buff, csv_file,
               TypedPrintConsumer(buff, field_types, fallback_type, other_type_handlers))
    seekstart(buff)
    read(buff, String)
end

export PrintConsumer


text:
using Vectorize
using Test
using Random
using Libdl
using SpecialFunctions
using LinearAlgebra

## Test Apple Accelerate
@static if Sys.isapple()
    include("AccelerateTests.jl")
else
    println("Accelerate not present. Aborting Accelerate tests")
end

## Test Yeppp
include("YepppTests.jl")

## Check for presence of VML and test VML
if Libdl.find_library(["libmkl_vml_avx"], ["/opt/intel/mkl/lib"]) != ""
    include("VMLTests.jl")
end

## Running tests over benchmarked Vectorized functions
include("VectorizeTests.jl")


text:
#=
Starting with the number 1 and moving to the right in a clockwise direction a 5 by 5 spiral is formed as follows:

21 22 23 24 25
20  7  8  9 10
19  6  1  2 11
18  5  4  3 12
17 16 15 14 13

It can be verified that the sum of the numbers on the diagonals is 101.

What is the sum of the numbers on the diagonals in a 1001 by 1001 spiral formed in the same way?
=#
#=
Notes:

#[1],2,[3],4,[5],6,[7],8,[9],10,11,12,[13],14,15,16,[17],18,19,20,[21],22,23,24,[25]
#1,3,5,7,9,13,17,21,25,31,37,43,49
#2, 4, 6
#3x3, 5x5, 7x7
=#
function calc()
  total,start = 1,1
  step = 2
  iterations = 3
  currentsquaresize = 3
  while currentsquaresize <= 1001
    start += step
    for i in start:step:(start+(iterations*step))
      total += i
    end
    start += (iterations*step)
    currentsquaresize += 2
    step += 2
  end
  total
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using TxKit

function run_test(registry, context, dim, image)
    new_method(registry, "debug") do method
        params = DebugParams(0.75)
        compute(context, method, image, Ref(params))

        download(image)

        map_read(image) do array
            display(array[1,:,:,1])
            println()

            display(array[1,:,:,2])
            println()

            display(array[1,:,:,3])
            println()
            @assert all(array[1,:,:,3] .== 0.0)

            display(array[1,:,:,4])
            println()
            @assert all(array[1,:,:,4] .== 0.75)
        end
    end
end

new_registry() do registry
    new_context(:gpu) do context
        dim = ImageDim(6, 3, 1, 4)

        new_image(:gpu, dim, Float32, 2, context) do image
            run_test(registry, context, dim, image)
        end
    end
end

# vim: ft=julia:ts=4:sw=4:et

max_stars_repo_name: vtavernier/txkit

max_issues_repo_name: vtavernier/txkit

max_forks_repo_name: vtavernier/txkit


text:
using Base.Test

import Base: ==

using KyotoCabinet

immutable K
  x::Int
end

immutable V
  a::Int
  b::String
end

==(x::K, y::K) = x.x == y.x
==(x::V, y::V) = (x.a == y.a) && (x.b == y.b)

KyotoCabinet.pack(v::ASCIIString) = convert(Array{Uint8,1}, v)
KyotoCabinet.unpack(T::Type{ASCIIString}, buf::Array{Uint8,1}) = bytestring(buf)

function KyotoCabinet.pack(k::K)
  io = IOBuffer()
  write(io, int32(k.x))
  takebuf_array(io)
end
function KyotoCabinet.unpack(T::Type{K}, buf::Array{Uint8,1})
  io = IOBuffer(buf)
  x = read(io, Int32)
  K(int(x))
end

function KyotoCabinet.pack(v::V)
  io = IOBuffer()
  write(io, int32(v.a))
  write(io, int32(length(v.b)))
  write(io, v.b)
  takebuf_array(io)
end
function KyotoCabinet.unpack(T::Type{V}, buf::Array{Uint8,1})
  io = IOBuffer(buf)
  a = read(io, Int32)
  l = read(io, Int32)
  b = bytestring(read(io, Uint8, l))
  V(int(a), b)
end

function test_get_set()
  file = tempdb()
  open(Db{K, V}(), file, "w+") do db
    db[K(1)] = V(1, "a")
    db[K(1999999999)] = V(2, repeat("b",100))
  end
  open(Db{K, V}(), file, "r") do db
    @assert V(1, "a") == db[K(1)]
    @assert V(2, repeat("b",100)) == db[K(1999999999)]
  end
end

function test_iter()
  file = tempdb()
  open(Db{K, V}(), file, "w+") do db
    db[K(1)] = V(1, "a")
    db[K(1999999999)] = V(2, repeat("b",100))
  end
  open(Db{K, V}(), file, "r") do db
    s0 = start(db)
    kv, s0 = next(db, s0)
    @assert K(1) == kv[1]
    @assert V(1, "a") == kv[2]

    kv, s0 = next(db, s0)
    @assert K(1999999999) == kv[1]
    @assert V(2, repeat("b", 100)) == kv[2]
  end
end

tempdb() = tempname() * ".kch"

test_get_set()
test_iter()


text:
@testset "examples" begin
    dir = "../examples"
    files = [ #"csg.jl",
             "example001.jl",
             "piping.jl",
             "radiused_union.jl",
             "fea.jl",
             "2d.jl"]
    for file in files
        @testset "$file" begin
                include("$dir/$file")
        end
    end
end


text:
struct SARAH_prox_iterable{R<:Real,C<:RealOrComplex{R},Tx<:AbstractArray{C},Tf,Tg}
    F::Array{Tf}            # smooth term - array of f_i s
    g::Tg                   # nonsmooth part
    x0::Tx                  # initial point
    N::Int                  # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of ∇f_i
    μ::Maybe{Union{Array{R},R}}  # convexity moduli of the gradients
    γ::Maybe{R}             # stepsize (a single scalar)
    m::Maybe{Int}           # number of inner loop updates
    ꞵ::Maybe{R}             # a scalar momentum
end

mutable struct SARAH_prox_state{R<:Real,Tx}
    γ::R                    # stepsize
    m::Maybe{Int}           # number of inner loop updates
    ꞵ::R                    # momentum
    av::Tx                  # average of ∇f_i s in outer loop
    w::Tx                   # w_{t-1} in the inner loop
    w_plus::Tx              # w_{t} in the inner loop
    ind::Array{Int}         # running idx set (1 to N)

    # some extra placeholders
    ∇f_temp::Tx             # placeholder for gradients
    temp::Tx                # placeholder for gradients
end

function SARAH_prox_state(γ::R, m, ꞵ::R, av::Tx, w::Tx, w_plus::Tx, ind) where {R,Tx}
    return SARAH_prox_state{R,Tx}(γ, m, ꞵ, av, w, w_plus, ind, copy(av), copy(av))
end

function Base.iterate(iter::SARAH_prox_iterable{R}) where {R}
    N = iter.N
    m = iter.m
    batch = 1
    ω = (3*(N-batch))/(2*batch*(N-1))
    ind = collect(1:N) # a vector of 1 to N
    # updating the stepsize
    if iter.γ === nothing
        γ = (2*sqrt(ω * m)) / (4*sqrt(ω * m) + 1)
    else
        γ = iter.γ # provided γ
    end

    if iter.ꞵ === nothing
        if iter.L === nothing
            @warn "smoothness parameter absent"
            return nothing
        else
            L_M = maximum(iter.L)
            ꞵ = 1/(L_M * sqrt(ω * m))
        end
    else
        ꞵ = iter.ꞵ
    end

    # initializing the vectors
    av = zero(iter.x0)
    w_plus = zero(iter.x0)
    for i = 1:N
        ∇f, ~ = gradient(iter.F[i], iter.x0)
        ∇f ./= N
        av .+= ∇f
    end
    w = copy(iter.x0)
    CIAOAlgorithms.prox!(w_plus, iter.g, w - γ * av, γ)
    w_plus .*= ꞵ
    w_plus .+= (1 - ꞵ) * w

    state = SARAH_prox_state(γ, m, ꞵ, av, w, w_plus, ind)
    return state, state
end

function Base.iterate(iter::SARAH_prox_iterable{R}, state::SARAH_prox_state{R}) where {R}
    # The inner cycle
    for i in rand(state.ind, state.m)
        gradient!(state.temp, iter.F[i], state.w_plus)
        gradient!(state.∇f_temp, iter.F[i], state.w)
        state.av .+= state.temp
        state.av .-= state.∇f_temp

        state.temp .= state.av
        state.temp .*= state.γ
        state.temp .*= -1
        state.temp .+= state.w_plus
        state.w .= state.w_plus
        CIAOAlgorithms.prox!(state.w_plus, iter.g, state.temp, state.γ)
        state.w_plus .*= state.ꞵ
        state.temp .= state.w
        state.temp .*= 1-state.ꞵ
        state.w_plus .+= state.temp
    end
    # full update
    state.w .= state.w_plus
    state.av .= zero(state.w)  # for next iterate
    for i = 1:iter.N
        gradient!(state.∇f_temp, iter.F[i], state.w)
        state.∇f_temp ./= iter.N
        state.av .+= state.∇f_temp
    end
    state.temp .= state.av
    state.temp .*= state.γ
    state.temp .*= -1
    state.temp .+= state.w
    CIAOAlgorithms.prox!(state.w_plus, iter.g, state.temp, state.γ)
    state.w_plus .*= state.ꞵ
    state.temp .= state.w
    state.temp .*= 1-state.ꞵ
    state.w_plus .+= state.temp

    return state, state
end

solution(state::SARAH_prox_state) = state.w_plus


text:
testfile = ARGS[1]

totab4(str) = " " ^ (div(length(str), 2) * 4)

function transform_line(line::AbstractString)
    m = match(r"^(.+)::String\s*,(.*)\n$", line)
    if m != nothing
        return transform_line("$(m.captures[1])::AbstractString,$(m.captures[2])\n")
    end

    m = match(r"^(.+)::String\s*\)(.*)\n$", line)
    if m != nothing
        return transform_line("$(m.captures[1])::AbstractString)$(m.captures[2])\n")
    end

    m = match(r"^(.+)::String\n$", line)
    if m != nothing
        return transform_line("$(m.captures[1])::AbstractString\n")
    end

    m = match(r"^(.+)<:\s+String\s*\}(.*)\n$", line)
    if m != nothing
        return transform_line("$(m.captures[1])<: AbstractString}$(m.captures[2])\n")
    end

    m = match(r"^(.+)<:\s+String\s*,(.*)\n$", line)
    if m != nothing
        return transform_line("$(m.captures[1])<: AbstractString,$(m.captures[2])\n")
    end

    m = match(r"^(.+)<:\s+String\s+$", line)
    if m != nothing
        return transform_line("$(m.captures[1])<: AbstractString\n")
    end

    m = match(r"^(.*)Dict\{\s*String\s*,(.*)\n", line)
    if m != nothing
        return transform_line("$(m.captures[1])Dict{AbstractString,$(m.captures[2])\n")
    end

    m = match(r"^(\s*)(.+)\{\s*String\s*\}(.*)\n", line)
    if m != nothing
        return transform_line("$(m.captures[1])$(m.captures[2]){AbstractString}$(m.captures[3])\n")
    end

    m = match(r"^(.+)\s+String\[\](.*)\n$", line)
    if m != nothing
        return transform_line("$(m.captures[1]) AbstractString[]$(m.captures[2])\n")
    end

    line
end

newlines = open(testfile, "r") do fh
    map(transform_line, readlines(fh))
end

# Delete so only max one empty line at end
while length(newlines) >= 2 && strip(newlines[end]) == "" && strip(newlines[end-1]) == ""
    deleteat!(newlines, length(newlines))
end

altname = tempname()
open(altname, "w") do fh
    print(fh, join(newlines, ""))
end

cmd = `diff $(testfile) $(altname)`
try
    run(cmd)
catch _e
end

# Remove file unless we should write it
if length(ARGS) < 2
    rm(altname)
elseif length(ARGS) >= 2 && in(lowercase(ARGS[2]), ["o", "-o", "--overwrite", "overwrite"])
    mv(altname, testfile; remove_destination = true)
else length(ARGS) >= 2 && in(lowercase(ARGS[2]), ["s", "-s", "--save", "write"])
    println("Altered file saved to: ", altname)
end


text:
using MoodleQuiz
using MoodleQuiz: moodle_xml_type
using Test

@testset "MoodleQuiz.jl" begin

    q = Question("shortanswer", 1, "", [ "What time is it?" ],
                 ShortAnswer("13:49"))
    @test moodle_xml_type(q) == "shortanswer"

    q = Question("numerical", 1, "", [ "1 + 1 =" ],
                 Numerical(2))
    @test moodle_xml_type(q) == "numerical"

    q = Question("essay", 1, "", [ "What did you do today?" ] )
    @test typeof(q.answer) == Essay
    @test moodle_xml_type(q) == "essay"

    q = Question("multilang essay", 1, "",
                 [
                     en"What time is it?"
                     cz"Kolik je hodin?"
                 ])

    q = Question("true/false", 1, "",
                 [ "Julia is cool." ],
                 True())
    @test typeof(q.answer) == TrueFalse
    @test moodle_xml_type(q) == "truefalse"

    q = Question("yes/no", 1, "",
                 [ "Is this question nice?" ],
                 Yes())
    @test moodle_xml_type(q) == "multichoice"

    q = Question("multichoice", 1, "",
                 [ "Moodle is (multiple of):" ],
                 MultiChoice(
                     Correct("LMS"),
                     Correct("web application"),
                     Wrong("Wiki"),
                     Wrong("desktop application"),
                 ))
    moodle_xml(q)
    @test moodle_xml_type(q) == "multichoice"

    q = Question("singlechoice", 1, "",
                 [ "Moodle is (single of):" ],
                 SingleChoice(
                     Correct("LMS"),
                     Wrong("Wiki"),
                     Wrong("desktop application"),
                 ))
    @test moodle_xml_type(q) == "multichoice"

    q = Question("cloze", 1, "",
                 [
                    "Moodle is"
                    ShortAnswer("LMS <&>")
                      "Its first version was release in year"
                      Numerical(2002)
                     "Is this quiz in Moodle?"
                     Yes()
                 ])
    @test moodle_xml_type(q) == "cloze"

    # The escaping rules are documented here:
    # https://docs.moodle.org/310/en/Embedded_Answers_(Cloze)_question_type#Detailed_syntax_explanations
    # but either the documentation or our implementation is not
    # complete, because Moodle fails to import the question if the
    # characters are in different order :-(
    q = Question("escaped cloze", 1, "",
                 [
                     "Characters that must be escaped in cloze sub-question are (write them space separated)"
                     ShortAnswer("# ~ \" \\ /")
                 ])
    @test occursin("\\# \\~ \\\" \\\\ \\/", MoodleQuiz.to_html(q))

    moodle_xml(q)
    moodle_xml(MoodleQuiz.default_category)
    moodle_xml(Quiz())

    if isinteractive()
        using LightXML
        save_file(moodle_xml(MoodleQuiz.default_category), "questions.xml")
    end
end


text:
using Revise
using Flipping
##
Directory_path = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/";
Exp_type = "Stimulations";
Exp_name = "DRN_Nac_Sert_ChR2"
Mice_suffix = "DN";
##
Directory_path = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/";
Exp_type = "Stimulations";
Exp_name = "DRN_Opto_Flipping"
Mice_suffix = "SD"
##
pokes, streaks, DataIndex = create_exp_dataframes(Directory_path,Exp_type,Exp_name, Mice_suffix)
##
DataIndex = Flipping.find_behavior(Directory_path,Exp_type,Exp_name, Mice_suffix)
findall(DataIndex[:Session].== "SD2_190131a.csv")
pokes, streaks, DataIndex = create_exp_dataframes(DataIndex[443:447,:])
pokes, streaks, DataIndex = create_exp_dataframes(DataIndex[1:10,:])
##
fn = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/run_task_photo/raw_data/SD2_190131a.csv"
check = process_pokes("/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/run_task_photo/raw_data/SD2_190131a.csv")
check[:,:PokeIn]
println(DataIndex[445,:])
union(DataIndex[:,:Session])
union(pokes[:,:Session])
findall(DataIndex[:,:Bhv_Path].== fn)
occursin("ddr","r")
DataIndex[:,:Session]
##
function tipi(df)
    res = DataFrame(Colonna = [], Tipo = [])
    for x in names(df)
        push!(res,(x,eltype(df[!,x])))
    end
    return sort!(res,:Colonna)
end
loaded = tipi(pokes_data)
done = tipi(pokes)
println((loaded,done))
check = DataFrame(Colonna = done[!,:Colonna], Tipo = done[!,:Tipo] .== loaded[!,:Tipo])
print(check)
##
DataIndex[1,:Preprocessed_Path]
##
DataIndex = Flipping.find_behavior(Directory_path, Exp_type, Exp_name,Mice_suffix)
##
Directory = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/run_task_2/All_LM"
Exp_name = "Learning"
Mice_suffix = "LM"
pokes, streaks, DataIndex = create_exp_dataframes(Directory)
DataIndex = Flipping.find_behavior(Directory)
union(DataIndex[:Day])
verify = by(DataIndex,:Day) do dd
    DataFrame(which = [SVector{size(dd,1),String}(union(dd[:MouseID]))])
end
println(verify)