text:
using Plots; pyplot()

include("../OceanParameters.jl")

## Replicate Bottom Loss Curves

θ₁ = π/4
@show BL = OceanParameters.bottom_loss(θ₁)

l = @layout [a b; c d]
θ₁ = deg2rad.(range(0, 90, length = 100))

cₚ = [1550. 1600. 1800.]
Numθ₁ = length(θ₁)
Numcₚ = length(cₚ)
BL = zeros(AbstractFloat, (Numθ₁, Numcₚ))
for nθ₁ = 1:length(θ₁)
	for ncₚ = 1:length(cₚ)
		BL[nθ₁, ncₚ] = OceanParameters.bottom_loss(θ₁[nθ₁]; cₚ = cₚ[ncₚ])
	end
end

p1 = plot(rad2deg.(θ₁), BL,
	xaxis = 0:30:90,
	yaxis = ("Loss (dB)", 0:5:15),
	label = cₚ,
	legendtitle = "cₚ (m/s)",
	legend = :bottomright)

αₚ = [1. 0.5 0]
Numθ₁ = length(θ₁)
Numαₚ = length(αₚ)
BL = zeros(AbstractFloat, (Numθ₁, Numαₚ))
for nθ₁ = 1:length(θ₁)
	for nαₚ = 1:length(αₚ)
		BL[nθ₁, nαₚ] = OceanParameters.bottom_loss(θ₁[nθ₁]; αₚ = αₚ[nαₚ])
	end
end

p2 = plot(rad2deg.(θ₁), BL,
	xaxis = 0:30:90,
	yaxis = 0:5:15,
	label = αₚ,
	legendtitle = "αₚ (dB/λ)",
	legend = :bottomright)

ρ₂ = [1.5e3 2e3 2.5e3]
Numρ₂ = length(ρ₂)
BL = zeros(AbstractFloat, (Numθ₁, Numρ₂))
for nθ₁ = 1:length(θ₁)
	for nρ₂ = 1:length(ρ₂)
		BL[nθ₁, nρ₂] = OceanParameters.bottom_loss(θ₁[nθ₁]; ρ₂ = ρ₂[nρ₂])
	end
end

p3 = plot(rad2deg.(θ₁), BL,
	xaxis = ("Grazing angle θ₁ (deg)", 0:30:90),
	yaxis = ("Loss (dB)", 0:5:15),
	label = ρ₂,
	legendtitle = "ρ₂ (kg/m³)",
	legend = :bottomright)

cₛ = [6e2 4e2 2e2 0.0]
Numcₛ = length(cₛ)
BL = zeros(AbstractFloat, (Numθ₁, Numcₛ))
for nθ₁ = 1:length(θ₁)
	for ncₛ = 1:length(cₛ)
		BL[nθ₁, ncₛ] = OceanParameters.bottom_loss(θ₁[nθ₁]; cₛ = cₛ[ncₛ], αₚ = 0.)
	end
end

p4 = plot(rad2deg.(θ₁), BL,
	xaxis = ("Grazing angle θ₁ (deg)", 0:30:90),
	yaxis = 0:5:15,
	label = cₛ,
	legendtitle = "cₛ (m/s)",
	legend = :bottomright)

pt = plot(p1, p2, p3, p4,
	layout = l)

savefig(pt, "img/BottomLoss_Parameters.png")

## Bottom Loss Types
θ₁ = deg2rad.(range(0., 90., length = Int(1e3)))
BottomTypes = ["Clay" "Silt" "Sand" "Gravel" "Moraine" "Chalk" "Limestone" "Basalt"]
ρ₂ = [1.5e3 1.7e3 1.9e3 2e3 2.1e3 2.2e3 2.4e3 2.7e3]
cₚ = [1.5e3 1575 1650 1.8e3 1950 2.4e3 3e3 5250]
cₛ = [5e1 0 0 0 6e2 1e3 1.5e3 2.5e3]
αₚ = [0.2 1. 0.8 0.6 0.4 0.2 0.1 0.1]
αₛ = [1.0 1.5 2.5 1.5 1.0 0.5 0.2 0.2]

NumBot = length(BottomTypes)
Numθ = length(θ₁)
BL = zeros(AbstractFloat, (Numθ, NumBot))
for nBot = 1:NumBot
	for nθ = 1:Numθ
		BL[nθ, nBot] = OceanParameters.bottom_loss(θ₁[nθ]; cₚ = cₚ[nBot], αₚ = αₚ[nBot], ρ₂ = ρ₂[nBot], cₛ = cₛ[nBot], αₛ = αₛ[nBot])
	end
end

pt = plot(rad2deg.(θ₁), BL,
	title = "Bottom Loss Types",
	yaxis = ("Bottom Reflection Loss (dB)", 0:5:20),
	xaxis = ("Grazing Angle (deg)", 0:15:90),
	label = BottomTypes,
	ylims = (0, 20))

savefig(pt, "img/BottomLoss_Types.png")

##


text:
using VoxelwiseConstantMapping
using Statistics: mean
using LinearAlgebra: norm
using Test

include("vxlwise_const.jl")
include("integral.jl")

ext: jl


text:
using Test

function compare_optimizer_states(optimizer1::OptimizerCmaEsDeap, optimizer2::OptimizerCmaEs, tolerance)

    @test optimizer1.lambda_ ≈ optimizer2.lambda_ atol = tolerance
    @test optimizer1.dim ≈ optimizer2.dim atol = tolerance
    @test optimizer1.chiN ≈ optimizer2.chiN atol = tolerance
    @test optimizer1.mu ≈ optimizer2.mu atol = tolerance
    @test optimizer1.weights ≈ optimizer2.weights atol = tolerance
    @test optimizer1.mueff ≈ optimizer2.mueff atol = tolerance
    @test optimizer1.cc ≈ optimizer2.cc atol = tolerance
    @test optimizer1.cs ≈ optimizer2.cs atol = tolerance
    @test optimizer1.ps ≈ optimizer2.ps atol = tolerance
    @test optimizer1.pc ≈ optimizer2.pc atol = tolerance
    @test optimizer1.centroid ≈ optimizer2.centroid atol = tolerance
    @test optimizer1.update_count ≈ optimizer2.update_count atol = tolerance
    @test optimizer1.ccov1 ≈ optimizer2.ccov1 atol = tolerance
    @test optimizer1.ccovmu ≈ optimizer2.ccovmu atol = tolerance
    @test optimizer1.C ≈ optimizer2.C atol = tolerance
    @test optimizer1.sigma ≈ optimizer2.sigma atol = tolerance
    @test optimizer1.damps ≈ optimizer2.damps atol = tolerance
    @test optimizer1.diagD ≈ optimizer2.diagD atol = tolerance
    @test optimizer1.B ≈ optimizer2.B atol = tolerance
    @test optimizer1.BD ≈ optimizer2.BD atol = tolerance
    @test optimizer1.genomes ≈ optimizer2.genomes atol = tolerance

end


text:
include("task_func.jl")

function aroundBorder!(r::Robot, side::HorizonSide; save = false)
    k = 0
    kk = 0
    while isborder(r, side)
        goToSideWithShag!(r, perSide(side), 1)
        k+=1
    end
    goToSideWithShag!(r, side, 1)
    while isborder(r, reverseSide(perSide(side)))
        goToSideWithShag!(r, side, 1)
        kk+=1
    end
    for _ = 1:k
        goToSideWithShag!(r, reverseSide(perSide(side)), 1)
    end
    if save
        return kk
    end
end

function drawPir(r)
    listOfShags = []
    side = West
    while true
        push!(listOfShags,goToSide!(r, side, save = true))
        if isRect(r, side)
            aroundBorder!(r, side)
        else
            break
        end
    end

    topCount = 0
    k = 0
    while !isborder(r, Nord)
        move!(r, Nord)
        topCount+=1
    end
    while !isborder(r, Ost)
        move!(r, Ost)
    end
    while !isborder(r, West)
        move!(r, West)
        k+=1
    end
    k+=1
    for _ = 1:topCount
        move!(r, Sud)
    end

    botCount = 0
    while !isborder(r, Sud)
        move!(r, Sud)
        botCount+=1
        localK = 0
        if k>0
            localK = 1
            putmarker!(r)
        end
        side = Ost
        while true
            while !isborder(r, side) && localK<k
                move!(r, side)
                putmarker!(r)
                localK+=1
            end

            if localK<k
                if isRect(r, side)
                    localK += aroundBorder!(r, side, save = true)
                    putmarker!(r)
                    localK+=1
                else
                    break
                end
            else
                break
            end
        end
        k-=1
        side = West
        while true
            while !isborder(r, side)
                move!(r, side)
            end

            if isRect(r, side)
                aroundBorder!(r, side)
            else
                break
            end
        end
    end

    for _ = 1:botCount
        move!(r, Nord)
    end

    for el = listOfShags
        goToSideWithShag!(r, Ost, el[1])
        if isborder(r, Ost)
            aroundBorder!(r, Ost)
        end
    end
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
function fixfiscal!(data::NCCSData; fiscaldatestr::String = FISCAL_DATE_FORMAT)
  #convert fiscal to dates
  rename!(data.df, :datefiscal=>:datefiscalstr)
  data.df[:datefiscal] =
    (s->ismissing(s) ? missing : Date(s,fiscaldatestr)).(data.df[:datefiscalstr])

  #correct bad fiscal dates
  data.df[:datefiscal] = ((fisyr,datefiscal)->
    (!ismissing(datefiscal)) && (Year(datefiscal) != fisyr) ?
      Date(fisyr, Dates.value(Month(datefiscal)), 1) : datefiscal
      ).(data.df[:fisyr], data.df[:datefiscal])

end

#this merges in a field or set of fields into data
#assumes the data being merged is logged returns
function mergelreturnfields!(data::NCCSData, indf::DataFrame, fieldstomerge::Vector{Symbol};
  fieldnames::Dict = Dict(f=>f for f ∈ fieldstomerge), #allows for renaming fields
  fieldtypes::Dict = Dict(f=>Union{Missing, eltype(indf[f])} for f ∈ fieldstomerge),
  fixfiscal::Bool = true, sortportfolio::Bool = true,
  lagfields::Bool = false) #allows for different field types

  #some time intensive activities
  fixfiscal && fixfiscal!(data) #only need to do this once
  sortportfolio && sort!(data.df, [:ein, :datefiscal])
  indf = sort(indf, :date)#sort!(indf, :date)

  #build an index for performance
  local indfindex::Dict{Union{Date,Missing},MInt} =
    Dict{Union{Date,Missing},MInt}(indf[i,:date]=>i for i ∈ 1:(size(indf,1)))
  indfindex[missing] = missing #provide a missing value if missing
  local idxbegin::Union{Int,Missing}
  local idxend::Union{Int,Missing}

  N::Int = size(data.df, 1)

  #pre-allocate the space
  for f ∈ fieldstomerge
    data.df[fieldnames[f]] = Vector{fieldtypes[f]}(undef, N)
    lagfields && (data.df[Symbol(:L, fieldnames[f])] = Vector{fieldtypes[f]}(undef, N))
  end

  fixfiscal && (data.df[:monthscovered] = Vector{MInt}(undef, N))

  #do merge at the NP level
  for subdf in groupby(data.df,:ein)
    subN::Int = size(subdf,1)

    for i::Int ∈ 1:subN
      idxend =  haskey(indfindex, subdf[i,:datefiscal]) ? indfindex[subdf[i,:datefiscal]] : missing
      if !ismissing(idxend)
        #assume eom dates (fiscal date given as yyyy-mm)
        #uses 12 month lag if no prior date is available

        if (((i==1) || (ismissing(subdf[i-1,:datefiscal]))) ||
          (subdf[i-1,:fisyr] != subdf[i,:fisyr] - 1))
            idxbegin = (idxend ≥ 12) ? idxend - 11 : missing
        elseif haskey(indfindex, subdf[i-1,:datefiscal])
            idxbegin = indfindex[subdf[i-1,:datefiscal]] + 1
        else
           idxbegin = missing
        end

        fixfiscal && (subdf[i, :monthscovered] = idxend - idxbegin + 1) #capture this, but only once

        #merge the data if it exists
        if !ismissing(idxbegin)
          for f ∈ fieldstomerge #merge the fields and the lag field if desired
            subdf[i,fieldnames[f]] = sum(indf[idxbegin:idxend, f])
            lagfields && (i≥2) && (subdf[i,Symbol(:L, fieldnames[f])] = subdf[i-1,fieldnames[f]])
          end
        end

      end
    end
  end

  return nothing
end

function mergeSP500!(data::NCCSData;
    refreshwrds::Bool = REFRESH_WRDS,
    minReturnCutoff = MIN_RETURN_CUTOFF, maxReturnCutoff = MAX_RETURN_CUTOFF)

  #get the S&P500 data and index it
  local spdf = getWRDSSP500(refreshwrds = REFRESH_WRDS)

  #preallocate and make a lookup table for efficiency
  data.df[:sp500ret] = Vector{MFloat64}(missing, size(data.df,1))
  data.df[:Lsp500ret] = Vector{MFloat64}(missing, size(data.df,1))
  data.df[:t30ret] = Vector{MFloat64}(missing, size(data.df,1))
  data.df[:lsp500t30] = Vector{MFloat64}(missing, size(data.df,1))
  data.df[:excess] = Vector{MFloat64}(missing, size(data.df,1))

  mergelreturnfields!(data, spdf, [:lvwretd, :lt30ret];
    fieldnames = Dict(:lvwretd=>:lsp500ret, :lt30ret=>:lt30ret),
    fixfiscal = true, sortportfolio=true, lagfields=true)

  data.df[:sp500ret] .= (exp).(data.df[:lsp500ret]) .- 1
  data.df[:Lsp500ret] .= (exp).(data.df[:Llsp500ret]) .- 1

  data.df[:t30ret] .= (exp).(data.df[:lt30ret]) .- 1
  data.df[:excess] .= data.df[:return] .- data.df[:t30ret]
  data.df[:lsp500t30] .= data.df[:lsp500ret] .- data.df[:lt30ret]

  #NOTE: Move these to another method???
  data.df[:pcontributions] = (data.df[:contributions]) ./ data.df[:lagnetassets]
  data.df[:pexpenses] = (data.df[:totexp]) ./ data.df[:lagnetassets]

  #trim rediculous values
  for v ∈ [data.df[:excess], data.df[:pcontributions], data.df[:pexpenses]]
    v .= (f->((ismissing(f)) || (f>maxReturnCutoff) || (f<minReturnCutoff)) ? missing : f).(v)
  end

  return data
end

#loads the Fama French Factors from a file
function loadFF(ffnamesuffix::String; refreshoutsidedata::Bool = true,
  ffpath::String = FF_PATH, ffname::String = FF_NAME,
  ffdateformat::String = FF_DATE_FORMAT)

  local ffdf::DataFrame
  local fffile::String = "$ffpath\\$(ffname)$(ffnamesuffix)"

  #load from a jls file if possible and desired
  if refreshoutsidedata || (!isfile("$fffile.csv"))
    ffdf= CSV.read("$fffile.csv")
    rename!(ffdf, :date=>:olddate)
    ffdf[:date] = (d->Date("$d",ffdateformat)).(ffdf[:olddate])
    deletecols!(ffdf, :olddate)

    serialize("$fffile.jls", ffdf)

  else
    ffdf = deserialize("$fffile.jls")
  end

  return ffdf
end

function mergeFF!(data::NCCSData, FFmodelsuffix::String; refreshoutsidedata::Bool = true)::Nothing
  N::Int = size(data.df,1)

  ffdf::DataFrame = loadFF(FFmodelsuffix, refreshoutsidedata=refreshoutsidedata)

  #merge in the raw data (fields from 3 and 5 factor models)
  fieldstomerge::Vector{Symbol} = FF_SUFFIX_FIELDS[FFmodelsuffix]
  lfieldstomerge::Vector{Symbol} = (s->Symbol(:l,s,FFmodelsuffix)).(fieldstomerge)

  ffN::Int = size(ffdf, 1)

  for i ∈ 1:length(fieldstomerge)
    ffdf[lfieldstomerge[i]] = Vector{MFloat64}(undef, ffN)
    ffdf[lfieldstomerge[i]] .= (f->ismissing(f) ? missing : log(1+f/100.0)).(ffdf[fieldstomerge[i]])
  end

  mergelreturnfields!(data, ffdf, lfieldstomerge, fixfiscal=false, lagfields=true)

  #readjust the mkt benchmark
  lmktt30sym::Symbol = Symbol("lmkt", FFmodelsuffix,"t30")
  data.df[lmktt30sym] = Vector{MFloat64}(undef, N)
  data.df[lmktt30sym] .= data.df[Symbol("lMkt_RF", FFmodelsuffix)] .+
    data.df[Symbol("lRF", FFmodelsuffix)] .- data.df[:lt30ret] #change rf benchmark

  return nothing
end

function mergeoutsidedatasources!(data::NCCSData; refreshoutsidedata::Bool = false,
  FFmodelsuffixes::Vector{String} = ["3", "5"])

  mergeSP500!(data::NCCSData)
  (FFmodelsuffix->mergeFF!(data, FFmodelsuffix)).(FFmodelsuffixes)
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
include("common.jl")

Δt = 0.01
mmesh = Mesh(30, 30, 1/30)
bdnode = bcnode(mmesh)
bdnode = [bdnode; bdnode .+ mmesh.ndof]
nu = Variable(0.001) * ones(get_ngauss(mmesh))
nodes = fem_nodes(mmesh)
u = @. sin(2π * nodes[:,1])
v = @. sin(2π * nodes[:,2])
u0 = [u;v]
u0[bdnode] .= 0.0
us = solve_burgers(u0, 10, nu)

U = matread("fenics/fwd1.mat")["U"]
loss = sum((us - U)^2)
sess = Session(); init(sess)
# @info run(sess, loss)
BFGS!(sess, loss)


text:
using NormalHermiteSplines
using DoubleFloats
using Test

@testset "NormalHermiteSplines.jl" begin

include("1D.jl")
include("2D.jl")
include("3D.jl")

end


text:
include("../src/transducers.jl")

using Debug, Transducers
@debug function check_compose()
    fn(x) = x^2
  @bp
  transduce(dedupe(tmap(fn)), tpush!, Any[], [1, 1, 2, 2, 3, 3])
end


text:
#this example works on Windows
ccall( (:GetDoubleClickTime, "User32"), stdcall,
	Uint, (), )

ccall( (:clock, "libc"), Int32, ())

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("ImageCreator.jl")
include("dipole.jl")

const c = 299792458.
const mu0 = 4*pi*1e-7
const eps0 = 1/(mu0*c^2)

#cavity
L=5
l=4
h=3
losses= 0.998
order=10

const freq=.5e9

##dipole
#dipole position
X=1.2
Y=1
Z=1.3
tilt=pi/2-acos(sqrt(2/3));
azimut=pi/4
phase=0
#dipole moment
#total time averaged radiated power P= 1 W dipole moment => |p|=sqrt(12πcP/µOω⁴)
Pow=1
amplitude=sqrt(12*pi*c*Pow/(mu0*(2*pi*freq)^4))

#Images
POS=IC(L,l,h,X,Y,Z,tilt,azimut,phase,amplitude,order)

#using HDF5,JLD
#@load "POS.jld" POS

numberofimages=1+2*order+(2*order*(order+1)*(2*order+1))/3
POS=POS[1:numberofimages,:]

#observation points
dx=0.02
dy=0.02
dz=0.02
const  nx=int(L/dx)
const  xmax=L
const  ny=int(h/dy)
const  ymax=l
const  x=linspace(0,L,nx)
const  z=1
const  y=linspace(0,l,ny)

println("Computing the radiation...")
Z=zeros(nx,ny)
Et=zeros(nx,ny)
Bt=zeros(nx,ny)
for i=1:nx
  perc=round(i/nx*1000)/10
  println("$perc %")
  for j=1:ny
    r=[x[i],y[j],z] #position of the receiver
    E=zeros(Complex128,3)
    B=zeros(Complex128,3)
    for m=1:numberofimages
      p=vec(POS[m,1:3]) #image dipole moment
      R=vec(POS[m,4:6]) #image dipole position
      ord=POS[m,7] #order of the dipole
      Ed,Bd=Hertz_dipole (r, p*losses^ord, R, phase, freq) #fields computation
      E+=Ed
      B+=Bd
    end
    Z[i,j]=sqrt(sum(abs(E).^2))/sqrt(sum(abs(B).^2))*mu0 #real(E).^2#0.5*numpy.cross(E.T,conjugate(B.T)) #impedance
    Et[i,j]=sqrt(sum(abs(E).^2)) #total E-field
    Bt[i,j]=sqrt(sum(abs(B).^2)) #total B-field
  end
end

#Some figures
using PyPlot

figure(1) #total E-field
title("\$\\vert E\\vert\$ (V/m)")
pcolor(x,y,Et',cmap="jet")
#clim(0,20)
colorbar()
axis("scaled")
xlim(0,L)
ylim(0,l)
grid()
xlabel("\$x\$/m")
ylabel("\$y\$/m")
savefig("E_$order.png",bbox="tight")
clf()

figure(2) #total H-field
title("\$\\vert H\\vert\$ (A/m)")
pcolor(x,y,Bt'/mu0,cmap="jet")
#clim(0,.2)
colorbar()
axis("scaled")
xlim(0,L)
ylim(0,l)
grid()
xlabel("x/m")
ylabel("\$y\$/m")
savefig("H_$order.png",bbox="tight")
clf()

figure(3) #wave impedance
title("\$Z\$ (\$ \\Omega\$)")
pcolor(x,y,Z',cmap="jet")
clim(0,1000)
colorbar(ticks=[0,100,200,300,120*pi,500,60*pi^2,700,800,900,1000])
axis("scaled")
xlim(0,L)
ylim(0,l)
grid()
xlabel("x/m")
ylabel("\$y\$/m")
savefig("Z_$order.png",bbox="tight")
clf()

max_stars_repo_name: manuamador/Z

max_issues_repo_name: manuamador/Z

max_forks_repo_name: manuamador/Z


text:
n = rand(1:10)
Z = Companion(randn(n))

#Special properties
@test Matrix(inv(Z)) ≈ inv(Matrix(Z))

#Matvec product
b = randn(n)
@test Z*b ≈ Matrix(Z)*b

m = rand(1:10)
A = randn(m, n)
@test A*Z ≈ A*Matrix(Z)


text:
using Weave
using Test

@testset "Weave" begin
    @testset "Chunk options" begin
        @info("Test: Chunk options")
        include("chunk_options.jl")
    end

    @testset "Error handling " begin
        @info("Testing error handling")
        include("errors_test.jl")
    end

    @testset "Eval in module" begin
        include("sandbox_test.jl")
    end

    @testset "Conversions" begin
        @info("Test: Converting")
        include("convert_test.jl")
    end

    @testset "Formatters" begin
        @info("Testing formatters")
        include("formatter_test.jl")
        include("markdown_test.jl")
        @info("Testing figure formatters")
        include("figureformatter_test.jl")
    end

    @testset "Rich output" begin
        @info("Testing rich output")
        include("rich_output.jl")
    end

    @testset "Plots" begin
        @info("Test: Weaving with Plots.jl")
        include("plotsjl_test.jl")
    end

    #@testset "Caching" begin
    #    @info("Test: Caching")
    #    include("cache_test.jl")
    #end
end

#@info("Test: Chunk options with Gadfly")
#include("chunk_opts_gadfly.jl")

#@info("Test: Weaving with Plots.jl")
#include("plotsjl_test.jl")
#include("publish_test.jl")


text:
export flip, flip!, randflip!

"""
    FlipStyle

Abstract type for flip styles.
"""
abstract type FlipStyle end

"""
    flip(site)

Defines how to flip this type of `site`.
"""
function flip end

flip(s::BT) where {BT <: BinarySite} = s == up(BT) ? down(BT) : up(BT)

"""
    flip!(S, I...)

Flips given configuration `S` at index `I...`.
"""
function flip! end

flip!(S::AbstractArray{BT, N}, I...) where {BT <: BinarySite, N} = flip!(S, I)
flip!(S::AbstractArray{BT, N}, I::Tuple) where {BT <: BinarySite, N} = (@inbounds S[I...] = flip(S[I...]); S)

"""
    randflip!(config) -> (proposed_index, config)
    randflip!(::FlipStyle, config) -> (proposed_index, config)

Flip the lattice configuration randomly using given flip style, default flip style is
[`UniformFlip`](@ref){1}, which choose **one** site in the configuration uniformly and flip it.

One should always be able to use `config[proposed_index]` to get the current value of
this lattice configuration. Whether one can change the site, depends on whether the
configuration is stored in a mutable type.
"""
function randflip! end

randflip!(S::AbstractArray) = randflip!(FlipStyle(S), S)

"""
    UniformFlip{N}

Choose `N` sites in the configuration uniformly and flip it
"""
struct UniformFlip{N} <: FlipStyle end
UniformFlip(N::Int) = UniformFlip{N}()
FlipStyle(S::AbstractArray{BT, N}) where {BT <: BinarySite, N} = UniformFlip(1) # use uniform(1) as default

function randflip!(::UniformFlip{1}, S::AbstractArray{BT}) where {BT <: BinarySite}
    proposed_index = rand(1:length(S))
    flip!(S, proposed_index)
    proposed_index, S
end

# TODO:
# 1. minimum sampling length of index can be calculated
#    use this to generate a shorter sequence
# 2. specialized method for static arrays (their length can be found in compile time)
function randflip!(::UniformFlip{N}, S::AbstractArray{BT}) where {BT <: BinarySite, N}
    proposed = zeros(Int, N)

    count = 0
    while count != N
        i = rand(1:length(S))
        i in proposed || (count += 1; proposed[count] = i)
    end

    @inbounds for i in 1:N
        flip!(S, proposed[i])
    end

    proposed, S
end


text:
repeat_stimulus(x, n_rep) = repeat(x, inner=(n_rep,1))

function simulate(model, n_type, max_rep)
    # generate random stimuli grouped by frequency of presentation
    stimuli = [random_stimulus(model, n_type) for _ in 0:max_rep]
    # study list in which stimuli are repeated 1 through max_rep times depending on group
    study_list = mapreduce(i->repeat_stimulus(stimuli[i+1], i), vcat, 1:max_rep)
    # encode study list
    encode!(model, study_list)
    # compute echo intensity for each frequency group
    return map(s->compute_intensity(model, s[1,:]), stimuli)
end


text:
using LinearAlgebra
using Distances

function kernel_BTG(r, theta)
    @assert length(theta) == 2
    return theta[1]^(r^(theta[2]))
end

# function kernel_gauss(d, x, y, θ)
#     @assert d == size(x, 1) "input vector x should in dimension d"
#     @assert d == size(y, 1) "input vector y should in dimension d"
#     @assert d == size(θ, 1) "length scale vector should in dimension d"
#     return exp(-dot(θ, (x-y).^2))
# end

function KernelMat(x0, x, kernel, param)
    n0 = size(x0, 1)
    n  = size(x, 1)
    dim = size(x, 2)
    Phi = zeros(n0, n)
    for i in 1:n0, j in 1:n
        if dim == 1
            Phi[i,j] = kernel(norm(x0[i] - x[j]), param)
        else
            Phi[i,j] = kernel(norm(x0[i,:] - x[j,:]), param)
        end
    end
    return Phi
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module DebugAdapter

import JuliaInterpreter

include("packagedef.jl")

end

max_stars_repo_name: fonsp/julia-vscode

max_issues_repo_name: fonsp/julia-vscode

max_forks_repo_name: fonsp/julia-vscode


text:
module Utilities

using LinearAlgebra # For dot
using OrderedCollections # for OrderedDict in UniversalFallback

import MutableArithmetics
const MA = MutableArithmetics

using MathOptInterface
const MOI = MathOptInterface

const MOIU = MOI.Utilities # used in macro

const SVF = MOI.VariableIndex
const VVF = MOI.VectorOfVariables
const SAF{T} = MOI.ScalarAffineFunction{T}
const VAF{T} = MOI.VectorAffineFunction{T}
const SQF{T} = MOI.ScalarQuadraticFunction{T}
const VQF{T} = MOI.VectorQuadraticFunction{T}

const VI = MOI.VariableIndex
const CI{F,S} = MOI.ConstraintIndex{F,S}

function print_with_acronym(io::IO, s::AbstractString)
    return print(io, replace_acronym(s))
end

function replace_acronym(s::AbstractString)
    s = replace(s, "MathOptInterface.Utilities" => "MOIU")
    s = replace(s, "MathOptInterface.Bridges" => "MOIB")
    s = replace(s, "MathOptInterface.Test" => "MOIT")
    s = replace(s, "MathOptInterface" => "MOI")
    return s
end

"""
    _reverse_dict(src::AbstractDict)

Reverse dictionary `src` so that values of the new dictionary are keys of `src`
and vice-versa. Also the values of `src` are assumed to be unique.

    _reverse_dict(dest::AbstractDict, src::AbstractDict)

Reverse dictionary so that values of `src` are key of `dest` and vice-versa.
`dest` must be empty. Also the values of `src` are assumed to be unique.
"""
function _reverse_dict end

include("CleverDicts.jl")
include("DoubleDicts.jl")

include("functions.jl")
include("mutable_arithmetics.jl")
include("sets.jl")
include("constraints.jl")
include("copy.jl")
include("results.jl")
include("variables.jl")

include("objective_container.jl")
include("variables_container.jl")
include("free_variables.jl")
include("vector_of_constraints.jl")
include("struct_of_constraints.jl")
include("model.jl")
include("sparse_matrix.jl")
include("product_of_sets.jl")
include("matrix_of_constraints.jl")
include("parser.jl")
include("mockoptimizer.jl")
include("cachingoptimizer.jl")
include("universalfallback.jl")
include("print.jl")

include("lazy_iterators.jl")

include("precompile.jl")
_precompile_()

end # module


text:
# Routine for estimating Lasso coefficients using GLMNet

function getLassoCoeff(outputs::Array{Float64,2},inputs::Array{Float64,2})

	numDependentVar::Int64 = size(outputs,2)
	numIndependentVar::Int64 = size(inputs,2)

	if(size(outputs,1) != size(inputs,1))
		throw(ErrorException("Number of output samples != number of input samples"))
	end

	coeff = zeros(Float64,numIndependentVar,numDependentVar)

	for d = 1:numDependentVar
		path_cv = glmnetcv(inputs[:,2:numIndependentVar],outputs[:,d]);
		path_ind::Int64 = indmin(path_cv.meanloss)
		coeff[2:numIndependentVar,d] = path_cv.path.betas[:,path_ind]
	end

	tmpMat = outputs - inputs[:,2:numIndependentVar]*coeff[2:numIndependentVar,:]
	for d = 1:numDependentVar
		coeff[1,d] = mean(tmpMat[:,d])
	end

	return coeff
end

max_stars_repo_name: rohitkannan/DD-SAA

max_issues_repo_name: rohitkannan/DD-SAA

max_forks_repo_name: rohitkannan/DD-SAA


text:
export G, C, M_P, K_B, SIGMA_T, SIGMA_SB,
M_SUN, YEAR_TO_SEC, SIGMA_E, H_PLANCK, ERG_TO_KEV, HZ_TO_KEV, UV_LOW_KEV, UV_HIGH_KEV

# all in cgs
const G = 6.6743e-8
const C = 29979245800.
const M_P = 1.67262192369e-24
const K_B = 1.380649e-16
const H_PLANCK = 6.62607015e-27
const SIGMA_T = 6.6524587321000005e-25
const SIGMA_E = SIGMA_T / M_P
const SIGMA_SB = 5.6703744191844314e-05
const M_SUN = 1.988409870698051e+33
const YEAR_TO_SEC = 31557600.0
const ERG_TO_KEV = 6.24150907e+08
const HZ_TO_KEV = 4.1356677e-18
const UV_LOW_KEV = 0.00387
const UV_HIGH_KEV = 0.06


text:
"""
A shape which contains all particles,``x > x_0``, necessary to simulate a plane-wave scattering from an infinite medium, for a reciever at the focal point, for time ``t < D / c``, where ``c`` is the sound speed of the background medium, and ``D`` is some chosen focal distance. The plane-wave travels towards the positive direction of the ``x`` axis.

More precisely, if the focal point is at ``(x_f,y_f)`` then the interior of the shape
is defined as
``(y - y_f)^2 < (D + x_0)^2 - x_f^2 - 2(D + x_0 - x_f)x``  and ``x > min(x_0, x_f)``
where ``D`` is the focal distance.
"""
struct TimeOfFlightPlaneWaveToPoint{T <: AbstractFloat,Dim} <: Shape{T,Dim}
    focal_point::AbstractVector{T}
    focal_distance::T
    minimum_x::T
end

TimeOfFlightPlaneWaveToPoint(focal_point::AbstractVector{T},  focal_distance::T;  minimum_x::T = zero(T)) where T <:AbstractFloat = TimeOfFlightPlaneWaveToPoint{T,length(focal_point)}(focal_point, focal_distance, minimum_x)

name(shape::TimeOfFlightPlaneWaveToPoint) = "Time of flight from planar source to the focal point"

# NOTE: this needs to be redone
function volume(shape::TimeOfFlightPlaneWaveToPoint{T}) where T <: AbstractFloat

    D = shape.focal_distance
    xf = shape.focal_point[1]
    x_min = shape.minimum_x
    x_max = (x_min + xf + D) / T(2)

    a = D + x_min + xf
    b = D + x_min - xf
    return 2/3*(sqrt(b*(a - 2*x_min)^3) - sqrt(b*(a - 2*x_max)^3))
end

import Base.issubset
function issubset(sphere::Sphere{T,Dim}, shape::TimeOfFlightPlaneWaveToPoint{T,Dim}) where {T, Dim}
    return (origin(sphere)[1] > shape.minimum_x) && (norm(origin(sphere) - shape.focal_point) <= shape.focal_distance - (origin(sphere)[1] - shape.minimum_x)  - T(2) * outer_radius(sphere))
end

function bounding_box(shape::TimeOfFlightPlaneWaveToPoint{T,Dim}) where {T,Dim}
    D = shape.focal_distance
    xf = shape.focal_point[1]
    y_f = shape.focal_point[2]
    x_min = shape.minimum_x
    x_max = (x_min + xf + D) / T(2)

    centre = zeros(Dim)
    centre[1] = (x_min + x_max) / T(2)
    centre[2] = y_f

    dimensions = ones(Dim) * T(2) * sqrt((D + x_min)^2 - xf^2 - 2*(D + x_min - xf)*x_min)
    dimensions[1] = x_max - x_min

    return Box(centre, dimensions)
end

function boundary_functions(shape::TimeOfFlightPlaneWaveToPoint{T,2}) where T

    D = shape.focal_distance
    xf = shape.focal_point[1]
    y_f = shape.focal_point[2]
    x_min = shape.minimum_x
    x_max = (x_min + xf + D) / T(2)

    function x(τ)
        check_boundary_coord_range(τ)

        if τ <= 1//3
            return (1 - 3*τ)*x_min + 3*τ*x_max
        elseif τ <= 2//3
            return (3*τ - 1)*x_min + (2 - 3*τ)*x_max
        else
            return x_min
        end
    end

    function y(τ)
        check_boundary_coord_range(τ)

        if τ <= 1//3
            return y_f + sqrt((D + x_min)^2 - xf^2 - 2*(D + x_min - xf)*((1 - 3*τ)*x_min + 3*τ*x_max))
        elseif τ <= 2//3
            return y_f - sqrt((D + x_min)^2 - xf^2 - 2*(D + x_min - xf)*((3*τ - 1)*x_min + 3*(2//3 - τ)*x_max))
        else
            return y_f + (6*τ - 5)*sqrt((D + x_min)^2 - xf^2 - 2*(D + x_min - xf)*x_min)
        end
    end

    return x, y
end


text:
using Sundials, Test, ForwardDiff
using Sundials: N_Vector, N_Vector_S
using LinearAlgebra

function mycopy!(pp, arr::Matrix)
    nj = size(arr, 2)
    ps = unsafe_wrap(Array, pp, nj)
    for j in 1:nj
        arr[:, j] = Sundials.asarray(ps[j])
    end
    arr
end

function mycopy!(arr::Matrix, pp)
    nj = size(arr, 2)
    ps = unsafe_wrap(Array, pp, nj)
    for j in 1:nj
        Sundials.asarray(ps[j])[:] = arr[:, j]
    end
end

f!(dy, t, y, p) = (dy[:] = y .* p)

function srhs(t, y, ydot, ys, ysdot)
    n = length(y)
    np = 2
    dyt = similar(y)
    chunk = min(n, 8)
    c1 = ForwardDiff.JacobianConfig(nothing, dyt, dyt, ForwardDiff.Chunk{chunk}())
    c2 = ForwardDiff.JacobianConfig(nothing, dyt, p, ForwardDiff.Chunk{chunk}())
    jac = ForwardDiff.jacobian((dy, y) -> f!(dy, t, y, p), dyt, y, c1, Val{false}())
    #jac = ReverseDiff.jacobian!(t1, y)
    ysdot[:] = jac * ys

    jac = ForwardDiff.jacobian((dy, p) -> f!(dy, t, y, p), dyt, p, c2, Val{false}())
    #jac = ReverseDiff.jacobian!(t2, p)
    ysdot[:, 1:np] += jac
end

"""
    sens(f!, t0, y0, p, tout; reltol, abstol)
Compute the solution and sensivities to the parametrized ODE problem defined by `f!(ẏ, t, y, p)`, starting at t0, y0, p.
Return the solutions `y` at times tout (>t0) as well as the corresponding sensivity matrices `ys`.
y[i, t] is the solutions component i at timestep t.
ys[i, j, t] is the i-th component sensivity wrt the j-th parameter at timestep t.
ys[i, np+j, t] the i-th component sensivity wrt the j-th initial condition value.
"""
function sens(
    f!::Function,
    t0::Float64,
    y0::Vector{Float64},
    p::Vector{Float64},
    tout::Vector{Float64};
    reltol::Float64 = 1e-5,
    abstol::Float64 = 1e-5,
)
    n = length(y0)
    np = length(p)
    ys0 = zeros(n, np .+ n)
    ys0[:, np .+ (1:n)] = Matrix(1.0I, n, n)

    #t1 = ReverseDiff.JacobianTape((dy,y)->f!(dy,0,y,p), dyt, y0)
    #t2 = ReverseDiff.JacobianTape((dy,p)->f!(dy,0,y0,p), dyt, p)

    pbar = abs.(vcat(p, y0))
    y, ys = cvodes(f!, srhs, t0, y0, ys0, p, reltol, abstol, pbar, tout)
end

### internals

## data structure dealing with the sundials callbacks

struct CVSData
    f::Any  # f(t,y,dy)
    fs::Any # fs()
    p::Any
    jys::Any
    jdys::Any
end

CVSData(f, fs, p, n::Int, nS::Int) =
    CVSData(f, fs, p, Array{Float64}(undef, n, nS), Array{Float64}(undef, n, nS))

function cvrhsfn(t::Float64, y::N_Vector, dy::N_Vector, data::CVSData)
    data.f(convert(Vector, dy), t, convert(Vector, y), data.p)
    return Sundials.CV_SUCCESS
end

function cvsensrhsfn(
    ns::Cint,
    t::Float64,
    y::N_Vector,
    dy::N_Vector,
    ys::N_Vector_S,
    dys::N_Vector_S,
    data::CVSData,
    tmp1::N_Vector,
    tmp2::N_Vector,
)
    jys = data.jys
    jdys = data.jdys
    mycopy!(ys, data.jys)
    data.fs(t, convert(Vector, y), convert(Vector, dy), jys, jdys)
    mycopy!(jdys, dys)
    return Sundials.CV_SUCCESS
end

## cvodes wrapper

"Given the sensivity problem, return (y,ys) where
y[i,t] is the solutions i-th componnent for timestep t and
ys[i,j,t] is the sensivity of the i-th component wrt to the j-th paramater, where
the last parameter indices correspond to the initial conditions components."
function cvodes(f, fS, t0, y0, yS0, p, reltol, abstol, pbar, t::AbstractVector)
    N, Ns = size(yS0)
    y = zeros(N, length(t))
    ys = zeros(N, Ns, length(t))
    tret = [t0]
    yret = similar(y0)
    ysret = similar(yS0)
    yS0n = [Sundials.NVector(yS0[:, j]) for j in 1:Ns]
    yS0nv = [N_Vector(n) for n in yS0n]
    pyS0 = pointer(yS0nv)
    crhs = Sundials.@cfunction(
        cvrhsfn,
        Cint,
        (Sundials.realtype, N_Vector, N_Vector, Ref{CVSData})
    )
    csensrhs = Sundials.@cfunction(
        cvsensrhsfn,
        Cint,
        (
            Cint,
            Sundials.realtype,
            N_Vector,
            N_Vector,
            N_Vector_S,
            N_Vector_S,
            Ref{CVSData},
            N_Vector,
            N_Vector,
        )
    )

    ##

    mem_ptr = Sundials.CVodeCreate(Sundials.CV_ADAMS)
    #mem_ptr = Sundials.CVodeCreate(Sundials.CV_BDF)
    cvode_mem = Sundials.Handle(mem_ptr)
    Sundials.CVodeSetUserData(cvode_mem, CVSData(f, fS, p, size(yS0)...))

    Sundials.CVodeInit(cvode_mem, crhs, t0, convert(N_Vector, y0))
    Sundials.CVodeSStolerances(cvode_mem, reltol, abstol)

    Sundials.CVodeSensInit(cvode_mem, Ns, Sundials.CV_STAGGERED, csensrhs, pyS0)
    Sundials.CVodeSetSensParams(cvode_mem, C_NULL, pbar, C_NULL)
    Sundials.CVodeSensEEtolerances(cvode_mem)
    for i in 1:length(t)
        @info "here1"
        Sundials.CVode(cvode_mem, t[i], yret, tret, Sundials.CV_NORMAL)
        @info "here2"
        Sundials.CVodeGetSens(cvode_mem, tret, pyS0)
        @info "here3"
        mycopy!(pyS0, ysret)
        y[:, i] = yret
        ys[:, :, i] = ysret
    end
    empty!(cvode_mem)
    y, ys
end

t0 = 0.0
t = [1.0, 2.0]
y0 = [1.0, 2.0]
p = [3.0, 4.0]
y, ys = sens(f!, t0, y0, p, t)
@test isapprox(y[1, 1], 20.0856, rtol = 1e-3)
@test isapprox(ys[2, 2, 2], 11924.3, rtol = 1e-3) # todo: check if these are indeed the right results


text:
function unpackfrac(q)
    numerator(q), denominator(q)
end

function add_one_variable(poly, newparent)
    R = parent(poly)
    base = base_ring(R)

    t = last(gens(newparent))

    polybuilder = MPolyBuildCtx(newparent)
    for (e, c) in zip(exponent_vectors(poly), coefficients(poly))
        push_term!(polybuilder, c, [e..., 0])
    end

    return finish(polybuilder), t
end

function erase_last_variable(poly, newparent)
    # assuming poly is indepent of last variable

    R = parent(poly)
    base = base_ring(R)

    polybuilder = MPolyBuildCtx(newparent)
    for (e, c) in zip(exponent_vectors(poly), coefficients(poly))
        push_term!(polybuilder, c, e[1:end-1])
    end

    return finish(polybuilder)
end

# TODO : todo
function change_parent_ring(poly, newparent)
    original = parent(poly)
    originalnvars = nvars(original)

    parentnvars = nvars(newparent)

    if originalnvars == parentnvars

    elseif originalnvars + 1 == parentnvars
        return add_one_variable(poly, newparent)
    elseif originalnvars == parentnvars + 1
        return erase_last_variable(poly, newparent)
    else
        @warn "failed to coerce polynomial $poly from $original to $parent"
    end
end

function unknown2known(u)
    libSingular.julia(libSingular.cast_number_to_void(u.ptr))
end

function julia(x::Singular.n_FieldElem{Singular.FieldElemWrapper{Nemo.GaloisField,Nemo.gfp_elem}})
    libSingular.julia(libSingular.cast_number_to_void(x.ptr)).data
end

function julia(x::Singular.N_Field{Singular.FieldElemWrapper{U, T}}) where {U, T}
    libSingular.julia(x.ptr).data
end

function julia(x)
    x
end

function singular2aa(poly::Singular.spoly{T}; base=false, new_ring=false) where {T}
    nvariables = length(gens(parent(poly)))
    xstrings = ["x$i" for i in 1:nvariables]
    if base == false
        base = base_ring(poly)
    end
    if new_ring == false
        new_ring, = AbstractAlgebra.PolynomialRing(base, xstrings, ordering=:lex)
    end
    change_base_ring(base, poly, parent=new_ring)
end

function double_singular2aa(
           poly::Singular.spoly{Singular.n_RingElem{Singular.spoly{T}}};
           base=false, new_ring=false) where {T}
    outer_change = singular2aa(poly)

    basebase = base_ring(parent(unknown2known(collect(coefficients(poly))[1])))

    nvariables = length(gens(parent(outer_change)))
    ystrings = ["y$i" for i in 1:nvariables]
    new_ring, = AbstractAlgebra.PolynomialRing(basebase, ystrings, ordering=:lex)

    inner_change = map_coefficients(
                      c -> singular2aa(unknown2known(c), base=basebase, new_ring=base),
                      outer_change)
    inner_change
end

function aa2singular(poly::MPoly{T}; base=false, new_ring=false) where {T}
    nvariables = length(gens(parent(poly)))
    xstrings = ["x$i" for i in 1:nvariables]
    if base == false
        base = base_ring(poly)
    end
    if new_ring == false
        new_ring, = Singular.PolynomialRing(base, xstrings, ordering=:lex)
    end
    change_base_ring(base, poly, parent=new_ring)
end

function double_aa2singular(poly::MPoly{T}; base=false, new_ring=false) where {T}
    nvariables = length(gens(parent(poly)))
    ystrings = ["y$i" for i in 1:nvariables]
    xstrings = ["x$i" for i in 1:nvariables]
    basebase = base_ring(base_ring(parent(poly)))
    if base == false
        base, = Singular.PolynomialRing(basebase, xstrings, ordering=:lex)
    end
    if new_ring == false
        new_ring, = Singular.PolynomialRing(base, ystrings, ordering=:lex)
    end
    change_base_ring(base_ring(new_ring), poly, parent=new_ring)
end

function Nemo.degree(f::AbstractAlgebra.Generic.Frac{T}) where {T}
    return max(degree(denominator(f)), degree(numerator(f)))
end

function Nemo.isconstant(f::AbstractAlgebra.Generic.Frac{T}) where {T}
    return isconstant(denominator(f)) && isconstant(numerator(f))
end

function Nemo.degrees(f::AbstractAlgebra.Generic.Frac{T}) where {T}
    return degrees(denominator(f)) + degrees(numerator(f))
end

function Base.length(x::Frac{T}) where {T}
    return max(length(numerator(x)), length(denominator(x)))
end

function applytofrac(f, x::Frac{T}) where {T}
    return map(f, (numerator(x), denominator(x)))
end

function tosingular(F::AbstractAlgebra.Generic.Rationals)
    Singular.QQ
end

function tosingular(F::Nemo.GaloisField)
    Singular.N_ZpField(Int(AbstractAlgebra.characteristic(F)))
end

function tosingular(F::Nemo.FlintRationalField)
    Singular.QQ
end

### TODO
function tosingular(F::Nemo.GFPMPolyRing)
    strings = string.(gens(F))
    ground = base_ring(F)
    Singular.PolynomialRing(tosingular(ground), strings)[1]
end

function tosingular(F::AbstractAlgebra.Generic.MPolyRing)
    strings = string.(gens(F))
    ground = base_ring(F)
    Singular.PolynomialRing(tosingular(ground), strings)[1]
end

function tosingular(F)
    F
end

function toaa(::Singular.Rationals)
    AbstractAlgebra.Generic.QQ
end

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

__randx = 10^4

function Base.rand(::AbstractAlgebra.Rationals{BigInt})
    return AbstractAlgebra.QQ(rand(1:__randx), rand(1:__randx))
end

function Base.rand(::Nemo.FlintRationalField)
    return Nemo.QQ(rand(1:__randx), rand(1:__randx))
end

function Base.rand(::Singular.Rationals)
    return Singular.QQ(rand(1:__randx), rand(1:__randx))
end

function Base.rand(wrapped::Singular.N_Field{Singular.FieldElemWrapper{U, T}}) where {U, T}
    return Base.rand( julia(wrapped) )
end

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


text:
module ExampleFs

using ..SQLiteFs

function lookup(req::FuseReq, parent::FuseIno, name::String, entry::CStructAccess{FuseEntryParam})
    println("lookup(req, $parent, $name) called")
    name != "xxx" && return Base.UV_ENOENT
    entry.ino = 4711
    entry.attr_timeout = 0.0
    entry.entry_timeout = 1e19
    st = entry.attr
    st.ino = 4712
    st.mode = 0o777
    st.uid = 1003
    st.gid = 10000
    st.size = 4712
    return 0
end

function getattr(req::FuseReq, ino, fi::CStruct{FuseFileInfo}, st::CStructAccess{Cstat}, attr_to::Ref{Float64})
    println("getattr(req, $ino, ...) called")
    st.ino = ino
    st.mode = 0o777
    st.uid = 10003
    st.gid = 10000
    st.size = 42
    return 0
end

args = ["mountpoint", "-d", "-h", "-V", "-f" , "-s", "-o", "clone_fd", "-o", "max_idle_threads=5"]

run() = main_loop(args, @__MODULE__)

end # module ExampleFs


text:
function assign_scenario(s::BasePriceScenario, config::VanillaOptions, params::VanillaParameters)
    params.cpricebase
end

function assign_scenario(s::OptimalPriceScenario, config::VanillaOptions, params::VanillaParameters)
    cprice_ubound = fill(Inf, config.N);
    cprice_ubound[1] = params.cpricebase[1];
    # Warning: If parameters are changed, the next equation might make base case infeasible.
    # If so, reduce tnopol so that we don't run out of resources.
    for i in 2:config.N
        if i > config.tnopol
            cprice_ubound[i] = 1000.0;
        end
    end
    cprice_ubound
end

function assign_scenario(s::Scenario, config::VanillaOptions, params::VanillaParameters)
    error("$(s) is not a valid scenario for v2013R (Vanilla flavour)");
end


text:
module Binned
import Base.Test
import Distributions
import StatsBase
const D = Distributions

using Entropy

function entropy()
    srand(1234)
    x = rand(1:5,1000)
    E = Entropy.entropy(x)
    Base.Test.@test_approx_eq E 2.3202534717146475
    println("Test passed. E = $E")
end

function nsb_entropy(ntrials::Integer,nvars::Integer)
	srand(1234)
	#Poisson variables
	X = rand(D.Poisson(0.2),(ntrials,nvars))
	E_true  = nvars*D.entropy(D.Poisson(0.2))
	M = (maximum(X)+1).^[0:nvars-1]
	N = X*M
	S = StatsBase.counts(N,minimum(N):maximum(N))
	return E_true,Entropy.NSBEntropy.simplefunc(N)
end

end


text:
module KMeans

using StatsBase

function initialize(colors, fixed, n)
    y = sample(1:size(colors)[2], n)
    x = sample(1:size(colors)[3], n)

    means = zeros(size(colors)[1], n)

    for i in 1:n
        means[:, i] = colors[:, y[i], x[i]]
    end

    return [fixed means]
end

function assign_clusters!(colors, means, clusters)
    nmeans = size(means)[2]

    distances = zeros((nmeans, size(colors)[2], size(colors)[3]))

    for i in 1:nmeans
        mean = means[:, i]

        distances[i, :, :] = sqrt.(sum((colors .- mean) .^ 2; dims=1))
    end

    clusters[:] = first.(Tuple.(argmin(distances; dims=1)))
end

function update_means!(colors, means, clusters, nfixed)
    nmeans = size(means)[2]

    for i in nfixed+1:nmeans
        mask = clusters .== i
        cluster = colors[:, mask]
        if size(cluster)[2] > 0
            means[:, i] = mean(cluster, dims=2)
        end
    end
end

function kmeans(colors, fixed, nmeans, nsteps)
    # Number of fixed means that do not get updated
    nfixed = if ndims(fixed) == 1
        size(fixed)[1]
    else
        size(fixed)[2]
    end

    # Initialize the means
    means = initialize(colors, fixed, nmeans)

    # Holds a cluster index per each color
    clusters = zeros(Int, size(colors)[2], size(colors)[3])

    # Main assign-update loop
    for _ in 1:nsteps
        assign_clusters!(colors, means, clusters)
        update_means!(colors, means, clusters, nfixed)
    end

    return means
end

end  # module


text:
__precompile__() # this module is safe to precompile
module SunAsAStar

include("common.jl")
export get_obs_loc, get_lst, get_solar_ra_dec, get_solar_hour_angle, get_earth_sun_dist, get_solar_info

# For interface to JPL Horizons & caching to disk
include("horizons.jl")
export Horizons

# For computing differential extinction
include("diff_solar_extinction.jl")
export DifferentialExtinction
export calc_Δv_diff_extinction

# For computing effects of apparent solar rotation rate
include("solar_rotation.jl")
export SolarRotation
export calc_Δfwhm_solar_rotation
#export CalculateFWHMDifference_SolarRotation_from_obs, CalculateFWHMDifference_SolarRotation_from_long_lat_alt
#export CalculateFWHMDifference_SolarRotation_from_loc_Ecliptic, CalculateFWHMDifference_SolarRotation_from_loc_Equatorial

# Temporary workaround since need to pin CSV
using CSV
export CSV

# Functions to precompile
precompile(SolarRotation.calc_Δfwhm_solar_rotation, (Float64,) )
precompile(DifferentialExtinction.calc_Δv_diff_extinction, (Float64,) )

end


text:
using StorageMirrorServer
using HTTP

using Tar
using Test
using Dates
using Random
using Suppressor

if VERSION < v"1.4"
    error("These tests require Julia at least v1.4")
end

@testset "StorageServer" begin
    include("tst_utils.jl")
    include("tst_server_utils.jl")
    include("tst_mirror_tarball.jl")
end


text:
"""
This package contains some codes from IterativeSolvers.jl: Copyright (c) 2013--2016 The Julia Language.
For the copyright notice, please check the LICENSE file.
"""
mutable struct BiCGStabIterable{
    opT, solT, vecT, matT, collectT,
    realT <: Real,
    scalarT <: Number}

    A::opT
    l::Int

    x::solT
    r_shadow::solT
    rs::collectT
    us::collectT

    maxiter::Int
    niter::Int
    tol::realT
    residual::realT

    gamma::vecT
    omega::scalarT
    sigma::scalarT
    M::matT

    function BiCGStabIterable(
        A::opT, l::Int, x::solT,
        r_shadow::solT, rs::collectT, us::collectT,
        maxiter::Int, niter::Int, tol::realT, residual::realT,
        gamma::vecT, omega::scalarT, sigma::scalarT,
        M::matT) where {opT,solT,vecT,matT,collectT,realT<:Real,scalarT<:Number}

        new{opT,solT,vecT,matT,collectT,realT,scalarT}(
            A,l,x,r_shadow,rs,us,maxiter,niter,tol,residual,gamma,omega,sigma,M)
    end
end

function init_bicgstab(
    x, A, b, l::Int=2;
    abstol::Real=zero(real(eltype(b))),
    reltol::Real=sqrt(eps(real(eltype(b)))),
    maxiter::Int=2000)

    T = typeof(x)
    N = size(x)
    niter = 0

    # Large vectors.
    r_shadow = T(undef, N)
    rand!(r_shadow)
    rs = Array{T}(undef, l+1)
    us = Array{T}(undef, l+1)
    for idx in 1:(l+1)
        rs[idx] = T(undef, N)
        us[idx] = T(undef, N)
        fill!(us[idx], zero(eltype(x)))
    end

    residual = rs[1]

    # Compute the initial residual = b - A * x
    mul!(residual, A, x)
    sub!(residual, b, residual)
    niter += 1

    gamma = zeros(eltype(x), l)
    omega = sigma = one(eltype(x))

    nrm = norm(residual)

    # For the least-squares problem
    M = zeros(eltype(x), l+1, l+1)

    # Stopping condition based on absolute and relative tolerance.
    tolerance = max(reltol * nrm, abstol)

    BiCGStabIterable(A, l, x, r_shadow, rs, us,
        maxiter, niter, tolerance, nrm,
        gamma, omega, sigma, M
    )
end

@inline converged(it::BiCGStabIterable) = it.residual ≤ it.tol
@inline start(::BiCGStabIterable) = 0
@inline done(it::BiCGStabIterable, iteration::Int) = it.niter ≥ it.maxiter || converged(it)

function iterate(it::BiCGStabIterable, iteration::Int=start(it))
    if done(it, iteration) return nothing end

    T = eltype(it.x)
    L = 2:(it.l + 1)

    it.sigma = -it.omega * it.sigma

    ## BiCG part
    for j = 1:(it.l)
        rho = dot(it.r_shadow, it.rs[j])
        beta = rho / it.sigma

        # us[:, 1 : j] .= rs[:, 1 : j] - β * us[:, 1 : j]
        for idx = 1:j
            sub!(it.us[idx], it.rs[idx], it.us[idx])
        end

        # us[:, j + 1] = Pl \ (A * us[:, j])
        mul!(it.us[j+1], it.A, it.us[j])

        it.sigma = dot(it.r_shadow, it.us[j+1])
        alpha = rho / it.sigma

        for idx = 1:j
            sub_scmul!(it.rs[idx], alpha, it.us[idx+1])
        end

        # rs[:, j + 1] = Pl \ (A * rs[:, j])
        mul!(it.rs[j+1], it.A , it.rs[j])

        # x = x + α * us[:, 1]
        add_scmul!(it.x, alpha, it.us[1])
    end

    # Bookkeeping
    it.niter += 2 * it.l

    ## MR part
    # M = rs' * rs
    for cidx in CartesianIndices(it.M)
        i1 = cidx[1]
        i2 = cidx[2]
        it.M[i1,i2] = dot(it.rs[i1], it.rs[i2])
    end

    # γ = M[L, L] \ M[L, 1]
    F = lu!(view(it.M, L, L))
    ldiv!(it.gamma, F, view(it.M, L, 1))

    for idx in L
        sub_scmul!(it.us[1], it.gamma[idx-1], it.us[idx])
        add_scmul!(it.x, it.gamma[idx-1], it.rs[idx-1])
        sub_scmul!(it.rs[1], it.gamma[idx-1], it.rs[idx])
    end

    it.omega = it.gamma[it.l]
    it.residual = norm(it.rs[1])

    it.residual, iteration + 1
end

"""
l: number of GMRES steps
"""
function bicgstab!(
    x, A, b, l::Int=2;
    abstol::Real=zero(real(eltype(b))),
    reltol::Real=sqrt(eps(real(eltype(b)))),
    maxiter::Int=2000,
    verbose::Bool=false)

    history = zeros(0)

    iterable = init_bicgstab(
        x, A, b, l;
        abstol=abstol, reltol=reltol,
        maxiter=maxiter)

    for (iteration, item) = enumerate(iterable)
        push!(history, iterable.residual)
        verbose && @printf("%3d\t%1.2e\n", iteration, iterable.residual)
    end

    if !converged(iterable)
        push!(history, -1)
        verbose && println("BiCGStab: not converged...")
    end
    # println()

    return history
end


text:
using ReduxImGui
using ReduxImGui.Redux
using ReduxImGui.CImGui

function naive_sliders(store::AbstractStore, get_state=Redux.get_state)
    if ReduxImGui.SliderInt(store, s->get_state(s).slider_ints["basic_slider_int"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
        v = ReduxImGui.get_value(get_state(store).slider_ints["basic_slider_int"])
        @info "Current value: $v"
    end
    CImGui.SameLine()
    ShowHelpMarker("CTRL+click to input value.")

    if ReduxImGui.SliderFloat(store, s->get_state(s).slider_floats["basic_slider_float"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
        v = ReduxImGui.get_value(get_state(store).slider_floats["basic_slider_float"])
        @info "Current value: $v"
    end

    if ReduxImGui.SliderFloat(store, s->get_state(s).slider_floats["basic_slider_float_curve"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
        v = ReduxImGui.get_value(get_state(store).slider_floats["basic_slider_float_curve"])
        @info "Current value: $v"
    end

    if ReduxImGui.SliderAngle(store, s->get_state(s).slider_angles["basic_slider_angle"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
        v = ReduxImGui.get_value(get_state(store).slider_angles["basic_slider_angle"])
        @info "Current value: $v"
    end

    if ReduxImGui.SliderString(store, s->get_state(s).slider_strings["basic_slider_string"])
        @info "This triggers $(@__FILE__):$(@__LINE__)."
        s = ReduxImGui.get_string(get_state(store).slider_strings["basic_slider_string"])
        @info "Current string: $s"
    end
    CImGui.SameLine()
    ShowHelpMarker("Using the format string parameter to display a name instead of the underlying integer.")
end


text:
#
#--------------------------------------#
export interpMat
#--------------------------------------#
"""
    Compute the Lagrange interpolation matrix from xi to xo.

    interpMat(xₒ,xᵢ)
"""
function interpMat(xo,xi)

    no = length(xo)
    ni = length(xi)

    a = ones(1,ni)
    for i=1:ni
        for j=1:(i-1)  a[i]=a[i]*(xi[i]-xi[j]); end
        for j=(i+1):ni a[i]=a[i]*(xi[i]-xi[j]); end
    end
    a = 1 ./ a # Barycentric weights

    J = zeros(no,ni)
    s = ones(1,ni)
    t = ones(1,ni)
    for i=1:no
        x = xo[i]
        for j=2:ni
            s[j]      = s[j-1]    * (x-xi[j-1]   )
            t[ni+1-j] = t[ni+2-j] * (x-xi[ni+2-j])
        end
        J[i,:] = a .* s .* t
    end

return J
end
#--------------------------------------#
export semInterp
#--------------------------------------#
"""
    Get local coordinate from physical coordinates

    u: field to interpolate
    xp,yp: physical coordinate of point
    rxe,rye,sxe,sye: inverse jacobian
    xe,ye: physical coordinates in element
"""
function semInterp(u,xp,yp,rxe,rye,sxe,sye,zr,zs,xe,ye
                ,tol=1e-8,maxiter=50)

# start from closest GLL point
d2 = @. (xe-xp)^2 + (ye-yp)^2
imin = argmin(d2)
r = zr[imin[1]]
s = zs[imin[2]]

Js = []
Jr = []
i = 0
while true
    Jr = interpMat(r,zr)
    Js = interpMat(s,zs)

    if(i == maxiter) break; end

    x = ABu(Js,Jr,xe)
    y = ABu(Js,Jr,ye)
    dx = xp .- x
    dy = yp .- y

    if(sqrt(dx.^2+dy.^2)[1] < tol) break; end

    # invere Jacobian
    rx = ABu(Js,Jr,rxe)
    ry = ABu(Js,Jr,rye)
    sx = ABu(Js,Jr,sxe)
    sy = ABu(Js,Jr,sye)

    # Newton iteration
    dr = rx * dx + ry * dx # [dr] = [rx ry] * [dx]
    ds = sx * dy + sy * dy # [ds]   [sx sy]   [dy]

    # add Hessian components
    #dr +=
    #ds +=

    println("r=",r,", s=",s,", iter=",i)

    if(abs(r+dr[1])>1. || abs(s+ds[1])>1.)
        println("approximation went to shit")
        return NaN
    end

    # update
    r += dr[1]
    s += ds[1]
    i += 1

end

return ABu(Js,Jr,u)[1]
end
#--------------------------------------#


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

using Markdown
using InteractiveUtils

# ╔═╡ 7a0d07f0-578a-11eb-359d-afc6edf042a0
begin
	in_pluto = @isdefined PlutoRunner
	if in_pluto
		using PlutoUI
		md"# Plots COVID19"
	end
end

# ╔═╡ bf11aa12-5788-11eb-2c7c-c979681e529d
begin
	using CSV
	using Dates
	using Statistics
	using DataFrames
	using Plots
end

# ╔═╡ 1031312c-5795-11eb-1df4-a71ffae9ba5d
begin
	if false
		using PlotlyJS
		plotlyjs()
	else
		gr()
	end
end;

# ╔═╡ e95448ec-578e-11eb-116f-8bde62ca2b8c
function moving_avarage(a::Vector, n::Int = 5)
	@assert length(a) ÷ 8 > n "n is to large for the array"
	[mean(a[max(1, i-n):min(i+n, end)]) for i in 1:length(a)]
end;

# ╔═╡ 5b05ded2-5795-11eb-020a-09c4f2bbc078
md"## UCI por región"

# ╔═╡ ead8c482-5788-11eb-1146-1f5095e4df4d
begin
	df_uci = CSV.read("../data/Datos-COVID19/output/producto8/UCI.csv", DataFrame)
end

# ╔═╡ 2be67732-578c-11eb-34a8-a799951cc203
begin
	dates = Date.(names(df_uci[4:end]));
	plot([], [], labels=nothing) # Empty plot
	for i in 1:size(df_uci, 1)
		uci_per_1000_people = 100_000 * Array(df_uci[i, 4:end]) / df_uci[i, 3]
		plot!(
			dates,
			moving_avarage(uci_per_1000_people, 5),
			labels=df_uci[i, 1],
			linewidth=1.5
		)
	end
	plot!(
		title="UCI por 100mil habitantes",
		legend=:outerright,
	)
end

# ╔═╡ 00285342-5794-11eb-30fb-cf6505a2a89f
begin
	plot(
		dates,
		moving_avarage(sum(Array.(eachrow(df_uci[:, 4:end])))),
		label=nothing,
		title="UCI Chile"
	)
end

# ╔═╡ 05010c5e-5796-11eb-39ce-e588ae1760d6
md"## Plan paso a paso"

# ╔═╡ 094cb4fc-5796-11eb-0d05-8356575e4b22
begin
	df_cuarentenas = CSV.read(
		"../data/Datos-COVID19/output/producto74/paso_a_paso.csv",
		DataFrame
	)
end

# ╔═╡ df7f2b0e-5796-11eb-3f9f-e5bdb7b07152
df_cuarentenas[df_cuarentenas["region_residencia"] .== "Arica y Parinacota", :]

# ╔═╡ 9d241c78-5797-11eb-25d8-2dcde74dd2df
# TODO: tratar de hacer unir esta tabla con lo de los casos actvos

# ╔═╡ Cell order:
# ╟─7a0d07f0-578a-11eb-359d-afc6edf042a0
# ╠═bf11aa12-5788-11eb-2c7c-c979681e529d
# ╠═1031312c-5795-11eb-1df4-a71ffae9ba5d
# ╠═e95448ec-578e-11eb-116f-8bde62ca2b8c
# ╟─5b05ded2-5795-11eb-020a-09c4f2bbc078
# ╠═ead8c482-5788-11eb-1146-1f5095e4df4d
# ╠═2be67732-578c-11eb-34a8-a799951cc203
# ╠═00285342-5794-11eb-30fb-cf6505a2a89f
# ╟─05010c5e-5796-11eb-39ce-e588ae1760d6
# ╠═094cb4fc-5796-11eb-0d05-8356575e4b22
# ╠═df7f2b0e-5796-11eb-3f9f-e5bdb7b07152
# ╠═9d241c78-5797-11eb-25d8-2dcde74dd2df

max_stars_repo_name: benjavicente/covid-maps

max_issues_repo_name: benjavicente/covid-maps

max_forks_repo_name: benjavicente/covid-maps


text:
module MatrixOptInterface

using LinearAlgebra
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

# export MatrixOptimizer

@enum ConstraintSense EQUAL_TO GREATER_THAN LESS_THAN

@enum VariableType CONTINUOUS INTEGER BINARY

include("matrix_input.jl")
include("allocate_load.jl")

end


text:
## subsystem initialization

export Shutdown, ismultithreaded

ismultithreaded() = convert(Core.Bool, API.LLVMIsMultithreaded())

for subsystem in [:Core, :TransformUtils, :ScalarOpts, :ObjCARCOpts, :Vectorization,
                  :InstCombine, :IPO, :Instrumentation, :Analysis, :IPA, :CodeGen, :Target]
    jl_fname = Symbol(:Initialize, subsystem)
    api_fname = Symbol(:LLVM, jl_fname)
    @eval begin
        export $jl_fname
        $jl_fname(R::PassRegistry) = API.$api_fname(R)
    end
end

## back-end initialization

export backends

const libllvm_backends = [:AArch64, :AMDGPU, :ARC, :ARM, :AVR, :BPF, :Hexagon, :Lanai,
                          :MSP430, :Mips, :NVPTX, :PowerPC, :RISCV, :Sparc, :SystemZ,
                          :VE, :WebAssembly, :X86, :XCore]

function backends()
    filter(libllvm_backends) do backend
        library = Libdl.dlopen(libllvm[])
        initializer = "LLVMInitialize$(backend)Target"
        Libdl.dlsym_e(library, initializer) !== C_NULL
    end
end

# generate subsystem initialization routines for every back-end
for backend in libllvm_backends,
    component in [:Target, :AsmPrinter, :AsmParser, :Disassembler, :TargetInfo, :TargetMC]

    initializer = "LLVMInitialize$(backend)Target"
    supported = :(Libdl.dlsym_e(Libdl.dlopen(libllvm[]), $initializer) !== C_NULL)

    jl_fname = Symbol(:Initialize, backend, component)
    api_fname = Symbol(:LLVM, jl_fname)
    @eval begin
        export $jl_fname
        function $jl_fname(error_on_use=true)
            if $supported
                ccall(($(QuoteNode(api_fname)),libllvm[]), Cvoid, ())
            elseif error_on_use
                error($"The $backend back-end is not part of your LLVM library.")
            end
        end
    end
end

# same, for initializing subsystems for all back-ends at once,
# reimplementes InitializeAll*s from the C-API which is `static inline`
for component in [:TargetInfo, :Target, :TargetMC, :AsmPrinter, :AsmParser, :Disassembler]
    jl_fname = Symbol(:Initialize, :All, component, :s)
    exprs = Expr[]
    for backend in libllvm_backends
        fname = Symbol(:Initialize, backend, component)
        push!(exprs, :($fname(false)))
    end

    @eval begin
        export $jl_fname
        function $jl_fname()
            $(exprs...)
            nothing
        end
    end
end

# same, for the native back-end
for component in [:Target, :AsmPrinter, :AsmParser, :Disassembler]
    jl_fname = Symbol(:Initialize, :Native, component)
    api_fname = Symbol(:LLVM, jl_fname)
    @eval begin
        export $jl_fname
        $jl_fname() = convert(Core.Bool, API.$api_fname()) &&
                      throw(LLVMException($"Could not initialize native $component"))
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ de66ce40-1945-11ec-19bd-53bbe9ed64c0
begin
    using Pkg
    Pkg.activate(".")
    using Revise
    using Javis
    using Colors
    using Images
    using PlutoUI
    using LaTeXStrings
    using LightXML
    using VideoIO
    using CoordinateTransformations
    using Rotations
	using BenchmarkTools
    import Luxor
end

# ╔═╡ aa4c48aa-e4fd-4cf6-b731-de97e55c3540
im = convert.(RGB, Luxor.@imagematrix begin
    background("blue")
    sethue("white")
    circle(O, 10, :fill)
    box(O, 50, 50, :stroke)
end 59 59)

# ╔═╡ ab0e4338-64cc-493c-8f71-cd77ad437ffb
im2 = convert.(RGB, Luxor.@imagematrix begin
    background("blue")
    sethue("white")
    circle(O, 10, :fill)
    # rotate(π/3)
    box(O, 50, 50, :stroke)
end 151 151)

# ╔═╡ 01735a86-406a-46a8-a3c6-d3b767eaee74
Javis.crop(im2, 59, 59) == im

# ╔═╡ a2878ded-81f8-4e7d-91da-b73ace797dda
function ground(color)
    (args...) -> begin
        Luxor.background(color)
        sethue("white")
    end
end

# ╔═╡ 38286288-86ff-4cee-9c58-296d0ed3c120
begin
    function ground_double_color(c11, c12, c21, c22, framenumber)
        if framenumber <= 50
            background(c11)
            sethue(c12)
        else
            background(c21)
            sethue(c22)
        end
    end
    myvideo = Video(200, 200)
    Background(1:100, (x, y, z) -> ground_double_color("blue", "black", "black", "red", z))
    circ1 = Object(1:100, JCircle(Point(0, 20), 20, action = :fill, color = "red"))
    # circ2 = Object(1:50, JCircle(Point(-20, -20), 20, action = :fill, color = "red"))
    act!(circ1, Action(anim_rotate(1π)))

    function postprocess_frame(frame_image, idx, frames)
        g(x) = RGBA(x.r, x.g, x.b, 0.6)

        frame_image = g.(frame_image)
        if idx <= 50
            frame_image
        else
            frame_image[size(frame_image, 1):-1:1, size(frame_image, 2):-1:1]
        end
    end

    function postprocess_frames_flow(frames)
        [frames; reverse(frames)]
    end
    tempdir = "temp3$(join(rand('a':'z', 10)))"

    v = render(
        myvideo,
        pathname = "output.gif",
        # tempdirectory="temp",
        postprocess_frame = postprocess_frame,
        # postprocess_frames_flow = postprocess_frames_flow,
        # liveview=true
    )
end

# ╔═╡ 00849422-a792-46c4-9623-06b3ee0a3750
N0f8

# ╔═╡ 48f6d415-6721-49bf-9390-fbae186f1c84
@btime let
    myvideo = Video(200, 200)
    Background(1:100, ground("black"))
    circ = Object(JCircle(Point(0, 20), 20, action = :fill, color = "white"))
    act!(circ, Action(anim_rotate(2π)))
    # tempdir = "temp3$(join(rand('a':'z', 10)))"

    render(
        myvideo,
        pathname = "output.gif",
    )
end

# ╔═╡ Cell order:
# ╠═de66ce40-1945-11ec-19bd-53bbe9ed64c0
# ╠═aa4c48aa-e4fd-4cf6-b731-de97e55c3540
# ╠═ab0e4338-64cc-493c-8f71-cd77ad437ffb
# ╠═01735a86-406a-46a8-a3c6-d3b767eaee74
# ╠═a2878ded-81f8-4e7d-91da-b73ace797dda
# ╠═38286288-86ff-4cee-9c58-296d0ed3c120
# ╠═00849422-a792-46c4-9623-06b3ee0a3750
# ╠═48f6d415-6721-49bf-9390-fbae186f1c84


text:
using NEAT
# output requests values to sum from X1, X2,...

# Create 2 input nodes
x0 = NEAT.BiasNode()
x1 = NEAT.InputNode(1.)
x2 = NEAT.InputNode(2.)

# Create 2 hidden nodes
h1 = NEAT.GeneNode(1)
h2 = NEAT.GeneNode(2)

# create connect from the inputs to hidden 1 and wire them up
w1 = NEAT.GeneConnection(x1,0,0,0)
w2 = NEAT.GeneConnection(x2,0,0,1)
push!(h1.xs, w1); w1.out = h1
push!(h1.xs, w2); w2.out = h1

# create connect from the inputs to hidden 2 and wire them up
w3 = NEAT.GeneConnection(x1,0,0,2)
w4 = NEAT.GeneConnection(x2,0,0,3)
push!(h2.xs, w3); w3.out = h2
push!(h2.xs, w4); w4.out = h2

# create connection from each hidden node to each output node
w5 = NEAT.GeneConnection(h1,0,0,3) # h1 ->  o1
w6 = NEAT.GeneConnection(h2,0,0,4) # h2 ->  o1
w7 = NEAT.GeneConnection(h1,0,0,5) # h1 ->  o2
w8 = NEAT.GeneConnection(h2,0,0,6) # h2 ->  o2

o1 = NEAT.GeneNode(5,[w5,w6],0.)
o2 = NEAT.GeneNode(6,[w7,w8],0.)

o1.f()
o2.f()

g = NEAT.Genome(2,2,[o1,o2],[x1,x2],[h1,h2],[w1,w2,w3,w4,w5,w6,w7,w8])
NEAT.addNode(g,2)


text:
using .SQLite
using .DBInterface

using Dates

const SQLITE_TYPE_MAPPINGS = Dict{Union{Type, Symbol}, Symbol}( # Julia => SQLite
  Char => :CHARACTER,
  String => :VARCHAR,
  Integer => :INTEGER,
  Int => :INTEGER,
  Float64 => :FLOAT,
  DateTime => :DATETIME,
  Time => :TIME,
  Date => :DATE,
  Bool => :BOOLEAN,
  Dict => :JSON
)

function database_column_type(dbtype::Type{SQLite.DB}, d::Union{Type, Symbol}) :: Symbol
    return SQLITE_TYPE_MAPPINGS[d]
end

function clean_table_query(table::Table, dbtype::Type{SQLite.DB})
    return "DELETE FROM $(table.name)"
end

database_kind(c::Type{SQLite.DB}) = Relational
close!(db::SQLite.DB) = DBInterface.close!(db)

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping


text:
module FEBases

using Gridap

export FEBasis

import Gridap: CellBasis
import Gridap: gradient
import Gridap: symmetric_gradient
import Base: div
import Gridap: trace
import Gridap: curl
import Gridap: inner
import Gridap: outer
import Gridap: varinner
import Base: +, -, *
import Gridap: restrict
import Gridap: Triangulation

struct FEBasis{B<:CellBasis,T<:Triangulation}
  cellbasis::B
  trian::T
end

function FEBasis(fespace::FESpace)
  b = CellBasis(fespace)
  trian = Triangulation(fespace)
  FEBasis(b,trian)
end

for op in (:+,:-,:(gradient),:(symmetric_gradient),:(div),:(trace),:(curl))
  @eval begin
    function ($op)(a::FEBasis)
      FEBasis($op(a.cellbasis),a.trian)
    end
  end
end

for op in (:+, :-, :*, :outer)
  @eval begin

    function ($op)(a::FEBasis,b::CellField)
      FEBasis($op(a.cellbasis,cellnewaxis(b,dim=1)),a.trian)
    end

    function ($op)(a::CellField,b::FEBasis)
      FEBasis($op(cellnewaxis(a,dim=1),b.cellbasis),b.trian)
    end

    function ($op)(a::FEBasis,b::Number)
      FEBasis($op(a.cellbasis,b),a.trian)
    end

    function ($op)(a::Number,b::FEBasis)
      FEBasis($op(a,b.cellbasis),b.trian)
    end

  end
end

function inner(a::FEBasis,b::CellField)
  varinner(a.cellbasis,b)
end

function varinner(a::FEBasis,b::CellField)
  varinner(a.cellbasis,b)
end

function inner(a::FEBasis,f::Function)
  b = CellField(a.trian,f)
  inner(a,b)
end

function inner(a::FEBasis,b::FEBasis)
  varinner(a.cellbasis,b.cellbasis)
end

function varinner(a::FEBasis,b::FEBasis)
  varinner(a.cellbasis,b.cellbasis)
end

Triangulation(a::FEBasis) = a.trian

function CellBasis(
  trian::Triangulation{D,Z},
  fun::Function,
  b::FEBasis,
  u...) where {D,Z}
  _u = [_prepare_cellbasis(ui) for ui in u]
  basis = CellBasis(trian,fun,b.cellbasis,_u...)
  FEBasis(basis,trian)
end

_prepare_cellbasis(ui) = ui

_prepare_cellbasis(ui::FEBasis) = ui.cellbasis

function restrict(feb::FEBasis,trian::BoundaryTriangulation)
  cb = restrict(feb.cellbasis,trian)
  FEBasis(cb,trian)
end

function restrict(feb::FEBasis,trian::SkeletonTriangulation)
  sp = restrict(feb.cellbasis,trian)
  _new_sp(sp,trian)
end

function _new_sp(sp::SkeletonPair{Z,T,N},trian) where {Z,T,N}
  b1 = FEBasis(sp.cellfield1,trian)
  b2 = FEBasis(sp.cellfield2,trian)
  SkeletonPair{Z,T,N}(b1,b2)
end

end # module


text:
using Test, LimnoSES

include("config.jl")
include("initialise.jl")
include("lake_state.jl")
include("models.jl")
include("reporting.jl")
include("evolve.jl")


text:
using PyCall
@pyimport skfmm as fmm

function perform_fast_marching(D, start_points)
    """
        Implementation of the fast marching algorithm in 2Ds using the skfmm library
        D : 2D weight matrix, must be positive
        start_points : 2 x n array, start_points[:,i] is the ith starting point
    """
    D_temp = D
    D_temp[start_points[1,:] + (start_points[2,:] - 1).*size(D)[1]] = 0
    return fmm.distance(D_temp) + 1e-15
end


text:
function calc_I_P(cur_reactor::AbstractReactor)
  if cur_reactor.is_pulsed
    cur_I = _calc_I_P_pulsed(cur_reactor)
  else
    cur_I = _calc_I_P_steady(cur_reactor)
  end

  cur_I *= cur_reactor.T_bar

  cur_I
end

function _calc_I_P_pulsed(cur_reactor::AbstractReactor)
  cur_G_PU = G_PU(cur_reactor)

  cur_G_V = G_V(cur_reactor)

  cur_numerator = cur_G_PU * ( K_BS(cur_reactor) )

  cur_denominator = cur_G_PU * ( 1 - K_CD(cur_reactor) * cur_reactor.sigma_v )

  cur_numerator += cur_G_V * K_VT(cur_reactor)

  cur_denominator -= cur_G_V * K_VI(cur_reactor)

  cur_numerator += G_CS(cur_reactor) * K_CS(cur_reactor)

  cur_I = cur_numerator

  cur_I /= cur_denominator

  cur_I
end

function _calc_I_P_steady(cur_reactor::AbstractReactor)
  cur_numerator = ( K_BS(cur_reactor) )

  cur_denominator = ( 1 - K_CD(cur_reactor) * cur_reactor.sigma_v )

  cur_I = cur_numerator

  cur_I /= cur_denominator

  cur_I
end


text:
## ExoplanetsSysSim/examples/profile.jl
## (c) 2015 Eric B. Ford

#workspace()
if !isdefined(:ExoplanetsSysSim)
  using ExoplanetsSysSim
end

module EvalModel
  export profile_model_eval
  using ExoplanetsSysSim

  sim_param_closure = SimParam()
  cat_phys_try_closure = KeplerPhysicalCatalog([])
  cat_obs_try_closure = KeplerObsCatalog([])
  summary_stat_try_closure =  CatalogSummaryStatistics()
  summary_stat_ref_closure =  CatalogSummaryStatistics()

  include(joinpath(Pkg.dir("ExoplanetsSysSim"), "src", "eval_model_macro.jl"))
  @make_evaluate_model(sim_param_closure,cat_phys_try_closure,cat_obs_try_closure,summary_stat_try_closure,summary_stat_ref_closure)

  function profile_model_eval(n::Integer)
    global sim_param_closure = setup_sim_param_demo()
    global cat_phys_try_closure = generate_kepler_physical_catalog(sim_param_closure)
    global cat_obs_try_closure = observe_kepler_targets_single_obs(cat_phys_try_closure,sim_param_closure)
    global summary_stat_ref_closure = calc_summary_stats_obs_demo(cat_obs_try_closure,sim_param_closure)

    param_base = make_vector_of_sim_param(sim_param_closure)

    evaluate_model_demo_scalar_param_scalar_ret(param::Float64) = sum(evaluate_model([param;param_base[2:end]]))

    io = open(joinpath(Pkg.dir("ExoplanetsSysSim"), "examples", "profile_syssim.tmp"),"w")
    Profile.print(io) # ,cols=800)
    close(io)
    #rm(io)

    Profile.init(10000000,0.001)
    eta_min = 0.5*get_real(sim_param_closure,"eta_pl")
    eta_max = 2.0*get_real(sim_param_closure,"eta_pl")
    #gc(); if VERSION >= VersionNumber(0,4,0)     gc_enable(false);   else     gc_disable();   end
    tic()
    for (i,eta) in enumerate(linspace(eta_min,eta_max,n) )
      local dist
      if i==2
         Profile.clear()
         Profile.clear_malloc_data()
      end
      @profile dist = evaluate_model_demo_scalar_param_scalar_ret(eta)
      println("eta= ",eta," d= ", dist)
    end
    toc()
    #gc_enable()
    pwd()
    io = open(joinpath(Pkg.dir("ExoplanetsSysSim"), "examples", "profile_syssim.out"),"w")
    Profile.print(io) # ,cols=800)
    close(io)
  end

end  # module EvalModel

using EvalModel

profile_model_eval(2)


text:
path = @get_scratch!("acrobot")

# generate residual methods
nq = acrobot_impact.nq
nu = acrobot_impact.nu
nc = acrobot_impact.nc
nz = nq + nc + nc
nz_nominal = nq
nθ = nq + nq + nu + 1

# Declare variables
@variables z[1:nz]
@variables z_nominal[1:nz_nominal]
@variables θ[1:nθ]
@variables κ[1:1]

# Residual
r_impact = residual(acrobot_impact, z, θ, κ)
r_impact = Symbolics.simplify.(r_impact)
rz_impact = Symbolics.jacobian(r_impact, z, simplify=true)
rθ_impact = Symbolics.jacobian(r_impact, θ, simplify=true)

# Build function
r_acrobot_impact_func = build_function(r_impact, z, θ, κ)[2]
rz_acrobot_impact_func = build_function(rz_impact, z, θ)[2]
rθ_acrobot_impact_func = build_function(rθ_impact, z, θ)[2]

rz_acrobot_impact_array = similar(rz_impact, Float64)
rθ_acrobot_impact_array = similar(rθ_impact, Float64)

@save joinpath(path, "impact.jld2") r_acrobot_impact_func rz_acrobot_impact_func rθ_acrobot_impact_func rz_acrobot_impact_array rθ_acrobot_impact_array
@load joinpath(path, "impact.jld2") r_acrobot_impact_func rz_acrobot_impact_func rθ_acrobot_impact_func rz_acrobot_impact_array rθ_acrobot_impact_array

# Residual
r_nominal = residual(acrobot_nominal, z_nominal, θ, κ)
r_nominal = Symbolics.simplify.(r_nominal)
rz_nominal = Symbolics.jacobian(r_nominal, z_nominal, simplify=true)
rθ_nominal = Symbolics.jacobian(r_nominal, θ, simplify=true)

# Build function
r_acrobot_nominal_func = build_function(r_nominal, z_nominal, θ, κ)[2]
rz_acrobot_nominal_func = build_function(rz_nominal, z_nominal, θ)[2]
rθ_acrobot_nominal_func = build_function(rθ_nominal, z_nominal, θ)[2]

rz_acrobot_nominal_array = similar(rz_nominal, Float64)
rθ_acrobot_nominal_array = similar(rθ_nominal, Float64)

path = @get_scratch!("acrobot")
@save joinpath(path, "nominal.jld2") r_acrobot_nominal_func rz_acrobot_nominal_func rθ_acrobot_nominal_func rz_acrobot_nominal_array rθ_acrobot_nominal_array
@load joinpath(path, "nominal.jld2") r_acrobot_nominal_func rz_acrobot_nominal_func rθ_acrobot_nominal_func rz_acrobot_nominal_array rθ_acrobot_nominal_array


text:
function twodQuadratureRule(rule::Int)
#= """
  twodQuadratureRule - calculate numerical integration points for triangular
                       elements

  Author: Jeff Borggaard, Virginia Tech
          part of FEMfunctions

  Licensing:
     This code is distributed under the MIT license.

  Usage:
  ```julia
    r,s,w = twodQuadratureRule(rule)
  ```

  Arguments:
  - `r`: xi coordinate of Gauss points
  - `s`: eta coordinate of Gauss points
  - `w`: Gauss weights corresponding to (r,s)

  Output arguments:
  - `rule`: Number of Gauss points:
""" =#

  r = zeros(Float64,rule)
  s = zeros(Float64,rule)
  w = zeros(Float64,rule)

  if rule == 1
    # The trivial linear triangle case
    r[1] = 1.0/3.0
    s[1] = 1.0/3.0

    w[1] = 0.5;

  elseif rule == 3
    # The following points correspond to a 3 point rule (quadratics)
    r[1] = 2.0/3.0;         s[1] = 1.0/6.0
    r[2] = 1.0/6.0;         s[2] = 2.0/3.0
    r[3] = 1.0/6.0;         s[3] = 1.0/6.0

    w[1] = 1.0/6.0
    w[2] = w[1]
    w[3] = w[1]

  elseif rule == 7
    # The following points correspond to a 7 point rule,
    # see Dunavant, IJNME, v. 21, pp. 1129-1148, 1995.
    # or Braess, p. 95.   (integrates 5th degree)

    t1 = 1.0/3.0;  t2 = (6.0+sqrt(15.0))/21.0;   t3 = 4.0/7.0 - t2

    r[1] = t1;            s[1] = t1
    r[2] = t2;            s[2] = t2
    r[3] = 1.0-2.0*t2;    s[3] = t2
    r[4] = t2;            s[4] = r[3]
    r[5] = t3;            s[5] = t3
    r[6] = 1.0-2.0*t3;    s[6] = t3
    r[7] = t3;            s[7] = r[6]

    t1 = 9.0/80.0;  t2 = (155.0+sqrt(15.0))/2400.0;  t3 = 31.0/240.0 - t2

    w[1]  = t1
    w[2]  = t2
    w[3]  = t2
    w[4]  = t2
    w[5]  = t3
    w[6]  = t3
    w[7]  = t3

  elseif rule == 13   # same Dunavant ref.  (integrates 7th degree)
    r[1]  = 0.0651301029022;  s[1]  = r[1]
    r[2]  = 0.8697397941956;  s[2]  = r[1]
    r[3]  = r[1];             s[3]  = r[2]
    r[4]  = 0.3128654960049;  s[4]  = 0.0486903154253
    r[5]  = 0.6384441885698;  s[5]  = r[4]
    r[6]  = s[4];             s[6]  = r[5]
    r[7]  = r[5];             s[7]  = r[6]
    r[8]  = r[4];             s[8]  = r[5]
    r[9]  = r[6];             s[9]  = r[4]
    r[10] = 0.2603459660790;  s[10] = r[10]
    r[11] = 0.4793080678419;  s[11] = r[10]
    r[12] = r[10];            s[12] = r[11]
    r[13] = 0.3333333333333;  s[13] = r[13]

    w[1]  = 0.0533472356088
    w[2]  = w[1]
    w[3]  = w[1]
    w[4]  = 0.0771137608903
    w[5]  = w[4]
    w[6]  = w[4]
    w[7]  = w[4]
    w[8]  = w[4]
    w[9]  = w[4]
    w[10] = 0.1756152574332
    w[11] = w[10]
    w[12] = w[10]
    w[13] =-0.1495700444677

    w = w/2.0

  elseif rule == 19   # same Dunavant ref.  (integrates 9th degree)
    r[1]  = 0.036838412054736; s[1]  = 0.741198598784498
    r[2]  = 0.741198598784498; s[2]  = 0.221962989160766
    r[3]  = 0.221962989160766; s[3]  = 0.036838412054736
    r[4]  = 0.741198598784498; s[4]  = 0.036838412054736
    r[5]  = 0.221962989160766; s[5]  = 0.741198598784498
    r[6]  = 0.036838412054736; s[6]  = 0.221962989160766
    r[7]  = 0.044729513394453; s[7]  = 0.910540973211095
    r[8]  = 0.044729513394453; s[8]  = 0.044729513394453
    r[9]  = 0.910540973211095; s[9]  = 0.044729513394453
    r[10] = 0.188203535619033; s[10] = 0.623592928761935
    r[11] = 0.188203535619033; s[11] = 0.188203535619033
    r[12] = 0.623592928761935; s[12] = 0.188203535619033
    r[13] = 0.437089591492937; s[13] = 0.125820817014127
    r[14] = 0.437089591492937; s[14] = 0.437089591492937
    r[15] = 0.125820817014127; s[15] = 0.437089591492937
    r[16] = 0.489682519198738; s[16] = 0.020634961602525
    r[17] = 0.489682519198738; s[17] = 0.489682519198738
    r[18] = 0.020634961602525; s[18] = 0.489682519198738
    r[19] = 0.333333333333333; s[19] = 0.333333333333333

    w[1]  = 0.043283539377289
    w[2]  = w[1]
    w[3]  = w[1]
    w[4]  = w[1]
    w[5]  = w[1]
    w[6]  = w[1]
    w[7]  = 0.025577675658698
    w[8]  = w[7]
    w[9]  = w[7]
    w[10] = 0.079647738927210
    w[11] = w[10]
    w[12] = w[10]
    w[13] = 0.077827541004774
    w[14] = w[13]
    w[15] = w[14]
    w[16] = 0.031334700227139
    w[17] = w[16]
    w[18] = w[16]
    w[19] = 0.097135796282799

    w = w/2.0

  else
    error("No quadrature rules other than 1, 3, 7, 13 or 19 implemented")
  end

  return r, s, w

end


text:
import YAML
config = YAML.load(open("config.yaml"))

using DiskJockey.constants
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.image
using DiskJockey.gridding
using HDF5

import PyPlot.plt
using LaTeXStrings

dv = DataVis(config["data_file"], true)[5]

# Conjugation is necessary for the SMA and ALMA
visibilities.conj!(dv)

qq = get_qq(dv)

# Find the indices that sort these in increasing order
ind = sortperm(qq)

qq = qq[ind]
uu = dv.uu[ind]
vv = dv.vv[ind]

fid = h5open("grid_VV.hdf5", "r")
# Load the sampled visibilities from HDF5
reals = read(fid, "real")
imags = read(fid, "imag")
VV = reals + imags .* im # Complex visibility

close(fid)

npoints = size(VV, 2)

mean_vis = mean(VV, 2)

fig, ax = plt[:subplots](nrows=3, figsize=(10, 9), sharex=true)

ms=0.1

for i=1:npoints
  vis = VV[:,i]
  ax[1][:plot](qq, abs(vis./mean_vis), ".", ms=ms)
  ax[2][:plot](qq, real(vis./mean_vis), ".", ms=ms)
  ax[3][:plot](qq, imag(vis./mean_vis), ".", ms=ms)
end

ax[3][:set_xlabel](L"$qq\quad k\lambda$")

ax[1][:grid]()
ax[2][:grid]()
ax[3][:grid]()

ax[1][:set_ylabel]("Mod")
ax[2][:set_ylabel]("Real")
ax[3][:set_ylabel]("Imag")

plt[:savefig]("visibilities.png", dpi=600)

max_stars_repo_name: iancze/DiskEmulator

max_issues_repo_name: iancze/DiskEmulator

max_forks_repo_name: iancze/DiskEmulator


text:
using TaylorIntegration, Elliptic
using Test
using LinearAlgebra: norm

@testset "Testing `bigfloats.jl`" begin

    local _order = 90
    local _abstol = 1.0E-77

    function pendulum!(dx, x, p, t) #the simple pendulum ODE
        dx[1] = x[2]
        dx[2] = -sin(x[1])
        nothing
    end

    @testset "Test ODE integration with BigFloats: simple pendulum" begin
        q0 = [big"1.3", 0.0] #the initial condition as a Vector{BigFloat}
        # T is the pendulum's librational period == 4Elliptic.K(sin(q0[1]/2)^2)
        # we will evaluate the elliptic integral K using TaylorIntegration.jl:
        G(x, p, t) = 1/sqrt(1-((sin(big"1.3"/2))^2)*(sin(t)^2)) # K elliptic integral kernel
        tvk, xvk = taylorinteg(G, 0.0, 0.0, BigFloat(π)/2, _order, _abstol)
        @test eltype(tvk) == BigFloat
        @test eltype(xvk) == BigFloat
        T = 4xvk[end] # T = 4Elliptic.K(sin(q0[1]/2)^2)
        @test typeof(T) == BigFloat
        @test T ≈ 4Elliptic.K(sin(q0[1]/2)^2) atol=eps(2.0) rtol=0.0

        t0 = 0.0 #the initial time

        tv, xv = taylorinteg(pendulum!, q0, t0, T, _order, _abstol; maxsteps=1)
        @test eltype(tv) == BigFloat
        @test eltype(xv) == BigFloat
        @test length(tv) == 2
        @test length(xv[:,1]) == 2
        @test length(xv[:,2]) == 2

        #note that T is a BigFloat
        tv, xv = taylorinteg(pendulum!, q0, t0, T, _order, _abstol)
        @test length(tv) < 501
        @test length(xv[:,1]) < 501
        @test length(xv[:,2]) < 501
        #the line below implies that we've evaluated the pendulum's period
        #up to an accuracy comparable to eps(BigFloat) ~ 1e-77!!!
        @test norm(xv[end,:].-q0,Inf) < 100eps(BigFloat)
    end
end


text:
# wrappers for functionality provided by the Metal library toolkit

include("metal/math.jl")


text:
using Documenter, IonBase

makedocs(;
    modules = [IonBase],
    format = Documenter.HTML(prettyurls = !("local" in ARGS)),
    pages = [
        "Home" => "index.md",
    ],
    repo = "https://github.com/Roger-luo/IonBase.jl",
    sitename = "IonBase.jl",
)

deploydocs(; repo = "github.com/Roger-luo/IonBase.jl")


text:
#=
https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl/blob/master/src/ReinforcementLearningZoo/src/algorithms/cfr/deep_cfr.jl
=#
function CounterfactualRegret.train!(sol::DeepCFRSolver, N::Int; show_progress::Bool = false, cb = () -> ())
    cb = Flux.Optimise.runall(cb)
    initialize!(sol)
    h0 = initialhist(sol.game)
    prog = Progress(N; enabled=show_progress)
    for t in 1:N
        cb()
        for p in 1:2
            for _ in 1:sol.traversals
                sol.traverser(sol, h0, p, sol.T)
            end
            train_value!(sol, p)
        end
        next!(prog)
        sol.T += 1.0f0
    end
    train_policy!(sol)
end

function initialize!(nn, init=Flux.glorot_normal)
    for p in Flux.params(nn)
        p .= init(size(p)...)
    end
end

"""
reset advantage networks and empty memory buffers
"""
function initialize!(sol::DeepCFRSolver)
    initialize!.(sol.V)
    initialize!(sol.Π)
end

function train_value!(sol::DeepCFRSolver, p::Int)
    initialize!(sol.V[p])
    opt = deepcopy(sol.advantage_opt)
    train_net!(
        sol.gpu,
        sol.V[p],
        sol.Mv[p].I,
        sol.Mv[p].r,
        sol.Mv[p].t,
        sol.batch_size,
        sol.value_batches,
        opt
    )
end

function train_policy!(sol::DeepCFRSolver, batches::Int=sol.strategy_batches, opt=sol.strategy_opt)
    train_net!(
        sol.gpu,
        sol.Π,
        sol.Mπ.I,
        sol.Mπ.σ,
        sol.Mπ.t,
        sol.batch_size,
        sol.strategy_batches,
        opt
    )
end

function train_net!(
    ::Val{true},
    dest_net,
    x_data,
    y_data,
    w,
    batch_size,
    n_batches,
    opt)

    isempty(x_data) && return nothing
    src_net = dest_net |> gpu

    input_size = length(first(x_data))
    output_size = length(first(y_data))
    p = Flux.params(src_net)

    _X = Matrix{Float32}(undef, input_size, batch_size)
    _Y = Matrix{Float32}(undef, output_size, batch_size)
    sample_idxs = Vector{Int}(undef, batch_size)
    idxs = 1:length(w)

    X = _X |> gpu
    Y = _Y |> gpu
    W = Vector{Float32}(undef, batch_size) |> gpu

    for i in 1:n_batches
        rand!(sample_idxs, idxs)
        fillmat!(_X, x_data, sample_idxs)
        fillmat!(_Y, y_data, sample_idxs)
        copyto!(X, _X)
        copyto!(Y, _Y)
        copyto!(W, w[sample_idxs])

        Loss = NetLoss(src_net, X, Y, W)

        gs = gradient(Loss, p)

        Flux.update!(opt, p::Flux.Params, gs)
    end

    Flux.loadmodel!(dest_net, src_net)
    nothing
end

function train_net!(
    ::Val{false},
    net,
    x_data,
    y_data,
    w,
    batch_size,
    n_batches,
    opt)

    isempty(x_data) && return nothing
    input_size = length(first(x_data))
    output_size = length(first(y_data))

    X = Matrix{Float32}(undef, input_size, batch_size)
    Y = Matrix{Float32}(undef, output_size, batch_size)
    W = Vector{Float32}(undef, batch_size)
    sample_idxs = Vector{Int}(undef, batch_size)
    idxs = 1:length(w)

    Loss = NetLoss(net, X, Y, W)
    p = Flux.params(net)

    for i in 1:n_batches
        rand!(sample_idxs, idxs)
        fillmat!(X::Matrix{Float32}, x_data, sample_idxs)
        fillmat!(Y::Matrix{Float32}, y_data, sample_idxs)
        copyto!(W::Vector{Float32}, @view w[sample_idxs])

        gs = gradient(Loss, p)

        Flux.update!(opt, p::Flux.Params, gs)
    end

    nothing
end

train_net!(args...) = train_net!(Val(false), args...)

function fillmat!(mat::AbstractMatrix, vecvec::AbstractVector, idxs)
    @inbounds for i in 1:size(mat, 2)
        mat[:,i] .= vecvec[idxs[i]]
    end
    return mat
end

"""
Weighted mean squared error
"""
wmse(ŷ,y,w) = sum(abs2.(ŷ .- y)*w ./ length(w))

struct NetLoss{NN,M,WGT}
    net::NN
    X::M
    Y::M
    W::WGT
end

(n::NetLoss)() = wmse(n.net(n.X), n.Y, n.W)


text:
module Rank2BinaryForms

using Combinatorics:combinations
using DocStringExtensions:SIGNATURES

# Use at least Nemo 0.24.0
# Install with
# ] add https://github.com/Nemocas/Nemo.jl.git
using Nemo

export crossratio, affineGamma, Image_affineGamma, set_getting_rootsofunity, buildpolynomial

#################################
#################################
# Roots of unity
#################################
#################################

rootsofunity(d) = k -> Nemo.root_of_unity(QQBar, d, k)

const ROOTS_OF_UNITY = Ref{Function}(rootsofunity)

"""
$(SIGNATURES)

Sets the method used to compute the `d`-th roots of unity. The input is a function `f` computing the roots of unity accepting calls as `f(d)(k)` and returning the `k`-th root of the `d`-th roots of unity. The order of the `d`-th roots of unity can be clockwise or anti-clockwise, but they must be ordered.
Examples of valid inputs are the following functions `f` and `g`:

```julia
    g(d) = k -> exp(2*pi*im*k/d)
    using Nemo
    CC = ComplexField(64);
    f(d) = k -> Nemo.root_of_unity(CC, d)^k;
```
"""
function set_getting_rootsofunity(f)
    ROOTS_OF_UNITY[] = f
end

## Other ways to get the roots of unity
# CC = ComplexField(64)
# g(d) = k -> Nemo.root_of_unity(CC, d)^k
# set_getting_rootsofunity(g)
#
# f(d) = k -> exp(2*pi*im*k/d)
# set_getting_rootsofunity(f)

#################################
#################################
# Cross ratio
#################################
#################################
"""
$(SIGNATURES)

Returns the cross ratio of four numbers.
There is also the method:

    crossratio(v) = crossratio(v...)
"""
crossratio(A, B, C, D) = (A - B) * (C - D) * inv((A - D) * (C - B))
crossratio(v) = crossratio(v...)

#################################
#################################
# Triplets and 4-tuples
#################################
#################################

"""
$(SIGNATURES)

Returns a vector of triplets `[i,j,k]` of integers `1 <= i < j < k <= d`, up to the action of the `d`-th dihedral group.
"""
triplets(d) = map(pair -> pushfirst!(pair, 1), combinations(2:d, 2))

"""
$(SIGNATURES)

Returns a vector of vectors resulting to `pushfirst` each integer `1:d` in `triplet`.
"""
tuples(d) = triplet -> [[i, triplet...] for i in 1:d if !(i in triplet)]
# Add zero and one as roots
# tuples(d) = triplet -> [[i, triplet...] for i in 1:d if !(i == last(triplet))]
# Add all the roots, if type has an infinity value.
# tuples(d) = triplet -> [[i, triplet...] for i in 1:d]

#################################
#################################
# Define Gamma
#################################
#################################
# Computing the i-th root of the polynomial image of Gamma
# @noinlne not needed
ithroot(d) = tuple -> crossratio(ROOTS_OF_UNITY[](d).(tuple))

# function affineGamma(d, triplet)
#     out = -1 .* ithroot(d).(tuples(d)(triplet))
#     return  [zero(first(out)), one(first(out)), out...]
# end
# affineGamma(d) = triplet -> affineGamma(d, triplet)

"""
$(SIGNATURES)

Returns a function that computes the unknown roots of the polynomial `Γ(triplet)`.
"""
affineGamma(d) = triplet -> ithroot(d).(tuples(d)(triplet))

"""
$(SIGNATURES)

Returns the roots of all the polynomials in the image of `Γ`.
"""
Image_affineGamma(d) = affineGamma(d).(triplets(d))

#################################
#################################
# Building the polynomials
#################################
#################################
"""
$(SIGNATURES)

Generic function to build the polynomials in the image of `Γ` from the list of roots.
"""
function buildpolynomial(roots, vars)
    x = first(vars)
    y = last(vars)
    out = x * y * (x + y)
    for r in roots
        out *= x - r * y
    end
    return out
end

function buildpolynomial(roots, vars::Vector{Symbol})
    x = first(vars)
    y = last(vars)
    out = Expr(:call, :*, x, y, Expr(:call, :+, x, y))
    for r in roots
        out = Expr(:call, :*, out, Expr(:call, :-, x, Expr(:call, :*, r, y)))
    end
    return out
end

buildpolynomial(vars) = roots -> buildpolynomial(roots, vars)

end


text:
include("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/scripts_for_figures/common_functions.jl")
include("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/scripts_for_figures/FK_sampling_functions.jl")

using DataFrames, DataFramesMeta, RCall, Distributions, StatsBase
R"library(copula)"

function R_q(q, b, k_n, alpha)
    t1 = exp(- ((b/k_n)^alpha)*((1 + q^(1/alpha))^alpha - q))
    t2 =  ((q^(1/alpha))/(1+q^(1/alpha)))^(1-alpha)
    return t1*t2
end

function sample_xi_n(alpha, b, theta, k_n, n)
   u =  Uniform(0,1)
   g = Gamma(theta/alpha +n, 1/((b/k_n)^(alpha)))
   r_n = rand(u,1)[1]
   dz_n  = rand(g,1)[1]
   while r_n > R_q(dz_n, b, k_n,alpha)
       u =  Uniform(0,1)
       r_n = rand(u,1)[1]
       dz_n  = rand(g,1)[1]
       if r_n<= dz_n
        # println("true")
       end
         #return (b/k_n)*((dz_n)^(1/alpha))
   end
   return (b/k_n)*((dz_n)^(1/alpha))
end

function density_xi_n(x, alpha, b, theta, k_n, n)
  term1 = exp(-((b/k_n + x)^alpha))
  term2 = (b/k_n + x)^(alpha - 1)
  term3 = x^(theta + (n-1)*alpha)
  return term1*term2*term3
end

## test xi_n density

#x_t = collect(range(0,1.2*10, length=100))
#xi = Array{Float64}(undef, 10000)
#for i in 1:10000
#  xi[i] = sample_xi_n(0.75, 1, 0.75, 1,1)
#end

## TEST for xi_n

#histogram(xi,bins=250)
#x_t = collect(range(0,1.2*10^2, length=1000))
#plot!(x_t,2000*map(x ->density_xi_n(x,0.75, 1.0, 0.75, 1, 1),x_t))

#histogram(sim.^(1/0.25),bins=250, normalize=true)
#plot!(x_t,map(x ->density_xi_n(x,0.25, 1.0, 0.25, 1, 1),x_t)/20)
#plot!(x_t,2000*map(x ->density_xi_n(x,0.75, 1.0, 0.75, 1, 1),x_t))

function SB_NGG(alpha_, b, theta,N)
  u_n =Array{Float64}(undef, N-1)
  k_n = 1
  for i in (1:(N-1))
    if i == 1
      k_n  = 1
    else
      k_n =prod(vec(ones(1,i-1)).-u_n[1:(i-1)])
    end
      xi_n = sample_xi_n(alpha_, b, theta, k_n,i)
    #println(xi_n)
    R"zn = rgamma(1,1- $alpha_, $b/$k_n + $xi_n)"
    #g = Gamma(1- alpha_, 1/(b/k_n + xi_n))
    #z_n  = rand(g,1)[1]
    z_n  = @rget zn
    #h_n  = b/k_n + xi_n
    R" xn <- retstable($alpha_, 1, h =$b/$k_n + $xi_n,method = 'MH')"
    x_n = @rget xn
    u_n[i] = z_n/ (z_n + x_n)
 end

  p = Array{Float64}(undef, N)
  p[1] = u_n[1]
  for l in 2:(N-1)
    p[l] = prod(vec(ones(1,l - 1)).-u_n[1:(l - 1)])*u_n[l]
  end
  p[N]=prod(vec(ones(1,(N-1))).- u_n)
  return p
end

#n_tr=10
#p_k = SB_NGG(.25, 1,0.25,n_tr)
#plot(1:n_tr,p_k)

function Pkn_NGG_SB(n, β, σ, M; runs=10^4)
  theta = 1
  b = (σ*theta*β)^(1/σ)
  alpha_ = σ
    array_clust_num = Array{Int64}(undef, runs)
    for i in 1:runs
        weights_NGG_SB = SB_NGG(alpha_,b, theta, M)
        c = wsample(1:M, weights_NGG_SB, n, replace=true)
        n_c =  length(unique(c))
        array_clust_num[i] = n_c
    end
    return proportions(array_clust_num, n)
end

function Pkn_NGG_SB_1(n, β, σ, M; runs=10^4)
  theta = (σ*β)
  b = 1
  alpha_ = σ
    array_clust_num = Array{Int64}(undef, runs)
    for i in 1:runs
        weights_NGG_SB = SB_NGG(alpha_,b, theta, M)
        c = wsample(1:M, weights_NGG_SB, n, replace=true)
        n_c =  length(unique(c))
        array_clust_num[i] = n_c
    end
    return proportions(array_clust_num, n)
end

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
using Maker
using Base.Test

fname = "t.txt"
open(fname, "w") do f
    println(f, "hello")
end
@test isfile(fname)
Maker.rm(fname)
@test !isfile(fname)
Maker.rm(fname)
@test !isfile(fname)

Maker.file("t.txt")

Maker.file("t1.txt") do
    open("t1.txt", "w") do f
        println(f, "hello")
    end
end
Maker.file("t2.txt") do
    open("t2.txt", "w") do f
        println(f, "hello")
    end
end
Maker.task("create files", ["t1.txt", "t2.txt"])
Maker.clean(["t1.txt", "t2.txt"])

make("create files")
@test isfile("t1.txt")
@test isfile("t2.txt")

make("clean")
@test !isfile("t1.txt")
@test !isfile("t2.txt")


text:
using Hungarian
using LinearAlgebra
using Random
using Distributions

"""
Run the algorithm from "Maximum Quadratic Assignment Problem: Reduction
from Maximum Label Cover and LP-based Approximation Algorithm" by
Makarychev, Manokaran, and Sviridenko (MMS) to approximately solve the QAP.
"""
function qap_mms(A, B, obj)
	@assert obj == :max
	x, y = adams_johnson_linearization(A, B, :max, integer=false)
	perm = rounding_mms(A, B, x, y)
	P = Matrix(I, size(A)...)[perm,:]
	return P, perm
end

"""
LP rounding step of the MMS algorithm.
"""
function rounding_mms(A, B, x, y)
	N = size(A, 1)
	n2 = Int(floor(N/2))

	La = randperm(N)[1:n2]
	Lb = randperm(N)[1:n2]

	Ra = setdiff(1:N, La)
	Rb = setdiff(1:N, Lb)

	ϕ̃ = fill(-1, N)
	for u in La
		probs = x[u,Lb]
		push!(probs, 1 - sum(probs))
		probvec!(probs)
		p′ = rand(Categorical(probs))
		p = (p′ == length(probs)) ? -1 : p′
		ϕ̃[u] = p
	end

	ϕ = fill(-1, N)
	for p in Lb
		us = findall(==(p), ϕ̃)
		if isempty(us) continue end
		u = rand(us)
		ϕ[u] = p
	end

	La = findall(!=(-1), ϕ)
	if isempty(La)
		return rounding_mms(A, B, x, y)
	end

	W = [sum(A[u,v] * B[ϕ[u],q] for u in La) for v in Ra, q in Rb]
	W = maximum(W) .- W
	matching = Hungarian.munkres(W)
	ψ = fill(-1, N)
	for i in 1:Int(ceil(N/2))
		j = findfirst(matching[i,:] .== Hungarian.STAR)
		ψ[Ra[i]] = Rb[j]
	end

	perm = [(ϕ[i] > 0) ? ϕ[i] : (ψ[i] > 0) ? ψ[i] : -1 for i in 1:N]

	missing_ind = findall(==(-1), perm)
	missing_vals = shuffle(setdiff(1:N, perm))
	perm[missing_ind] = missing_vals

	return perm
end

function probvec!(xs::Array{<:Real,1})
	z = zero(xs[1])
	s = z
	for i in 1:length(xs)
		if xs[i] < z
			xs[i] = z
		end
		s += xs[i]
	end
	if s != 1
		for i in 1:length(xs)
			xs[i] = xs[i] / s
		end
	end
	return xs
end


text:
push!(LOAD_PATH, "/mnt/g/Python Programme/BALL.jl/src" )
using BALL
using BenchmarkTools

pdb_filepath = "/mnt/g/Python Programme/BALL.jl/1,2-Benzodiazepine.pdb"

comp = System(pdb_filepath)

print(@benchmark MolecularGraph(comp))


text:
using Test #, Random

@testset "Testing solution to Exercise 1" begin

@testset "Running ex1.jl" begin
   include("../ex1.jl")
end;

@testset "Testing that variables exist" begin
   @test @isdefined(response_1a)
   @test @isdefined(response_1b)
   @test @isdefined(response_1c)
   @test @isdefined(response_1d)
   @test @isdefined(response_1e)
   @test @isdefined(response_1f)
   @test @isdefined(response_1g)
   @test @isdefined(response_1h)
   @test @isdefined(response_1i)
   @test @isdefined(response_1j)
end;

@testset "Testing that variables are not missing" begin
   @test check_type_isa(:response_1a,response_1a,Markdown.MD)
   @test check_type_isa(:response_1b,response_1b,Markdown.MD)
   @test check_type_isa(:response_1c,response_1c,Markdown.MD)
   @test check_type_isa(:response_1d,response_1d,Markdown.MD)
   @test check_type_isa(:response_1e,response_1e,Markdown.MD)
   @test check_type_isa(:response_1f,response_1f,Markdown.MD)
   @test check_type_isa(:response_1g,response_1g,Markdown.MD)
   @test check_type_isa(:response_1h,response_1h,Markdown.MD)
   @test check_type_isa(:response_1i,response_1i,Markdown.MD)
   @test check_type_isa(:response_1j,response_1j,Markdown.MD)
end;

@testset "Testing gradient calculations" begin
   @test maximum(abs.(grad_banana_2d_at_origin.-[0.8, -6.4]))<1e-5
   @test maximum(abs.(grad_banana_2d_at_minimum.-zeros(2)))<1e-5
end;

end; # Exercise 1


text:
module Internal

using Pkg
using UUIDs
using GitHub
using JuliaFormatter
# using TestEnv
using Comonicon
using GarishPrint
using Distributed
using Configurations
# using Glob: FilenameMatch, ismatch
using Comonicon.Tools: prompt

@static if Sys.which("git") === nothing # no system git
    using Git: git
else
    git() = Sys.which("git")
end

julia() = joinpath(Sys.BINDIR, Base.julia_exename())

include("utils/utils.jl")
include("blueprints/blueprints.jl")
# commands
include("commands/release.jl")
include("commands/clone.jl")
include("commands/compat.jl")
include("commands/create.jl")
include("commands/package.jl")
include("commands/format.jl")
include("commands/doc.jl")
# include("commands/test.jl")

# fix redirect color issue
# @static if !hasmethod(get, Tuple{Base.PipeEndpoint, Symbol, Any})
#     Base.get(::Base.PipeEndpoint, key::Symbol, default) = key === :color ? Base.get_have_color() : default
# end
# Base.get(::Base.PipeEndpoint, key::Symbol, default) = key === :color ? Base.get_have_color() : default

end