text:
# A FrequencyAdapter adapts the frequencies with which a set of
# values/strategies/methods should be applied/tried in an optimization problem.
# It is based on the Adaptive Coordinate Frequencies scheme described in:
#
# T. Glasmachers and U. Dogan, "Accelerated Coordinate Descent with
# Adaptive Coordinate Frequencies", 2013.
#
# but generalized so that it can support more than the adaptation of only
# the coordinates in a Coordinate Descent scheme. The things that are being
# adapted are identified by integers in 1:n, with n being the main parameter.
type FrequencyAdapter
  n::Integer
  eta::Float64
  c::Float64
  pmin::Float64
  pmax::Float64

  p::Vector{Float64}    # Current pi values.
  psum::Float64         # Current psum.
  a::Vector{Float64}    # Current ai values.
  deltahat::Float64     # Running average of the progress values.

  block::Vector{Int64}  # Current block of selected methods.

  numupdates::Int64     # Number of times we have been updated.
  min_updates::Int64    # Number of updates until we start adapting frequencies.

  FrequencyAdapter(n, c = 0.2, pmin = 0.05, pmax = 20) = begin
    eta = 1/n
    p = ones(Float64, n)/n
    new(n, eta, c, pmin, pmax,
      p, 1.0,
      zeros(Float64, n),
      0.0,
      shuffle(collect(1:n)),
      0, n)
  end
end

# Return the index of the next method that should be used. The base
# FrequencyAdapter first creates a block of randomly shuffled methods that should
# be applied and then selected the next one from it. If there is a current
# block which is non-empty, use that, if not create a new block. However,
# the first block is always random shuffles of all the methods since we need
# to learn about their effectiveness.
function next(fa::FrequencyAdapter)
  if length(fa.block) < 1
    create_new_block!(fa)
  end
  #println("Taking last from block = ", fa.block)
  return pop!(fa.block)
end

# Create a new block.
function create_new_block!(fa::FrequencyAdapter)
  block = Int64[]
  #print("Creating new block, psum = $(fa.psum), a = ", fa.a, ", p = ", fa.p)
  for(i in 1:fa.n)
    fa.a[i] += (fa.n * fa.p[i] / fa.psum)
    num_ai = convert(Int64, floor(fa.a[i]))
    if num_ai > 0
      block = vcat(block, ones(Int64, num_ai)*i)
      fa.a[i] -= num_ai
    end
  end
  # Due to rounding errors the block is sometimes empty so we select the one
  # with largest a value.
  if length(block) < 1
    #println("Empty block created. Rectifying., std(p) = ", std(fa.p))
    #print("psum = $(fa.psum), a = ", fa.a, ", p = ", fa.p)
    i = indmax(fa.a)
    fa.a[i] = 0
    block = [i]
  end
  #println("Created new block = ", block)
  fa.block = shuffle(block)
end

# Update the internal model of progress and success rate of each method based
# on the latest progress value of one method. Progress values should be larger
# the larger the progress/improvement was.
function update!(fa::FrequencyAdapter, methodIndex, progress)
  # If we already have collected a few samples of progress rates we can update
  # the pi. This is the common case.
  if fa.numupdates >= fa.min_updates
    unconstrained_pnew = fa.p[methodIndex] * exp(fa.c * (progress / fa.deltahat - 1))
    pnew = min(fa.pmax, max(fa.pmin, unconstrained_pnew))
    #print("i = ", methodIndex, ", pi = ", fa.p[methodIndex], ", pnew = ", pnew, ", psum = ", fa.psum)
    fa.psum += (pnew - fa.p[methodIndex])
    fa.p[methodIndex] = pnew
    #print(", new psum = ", fa.psum, ", deltahat = ", fa.deltahat)
    fa.deltahat = (1 - fa.eta) * fa.deltahat + fa.eta * progress
    #println(", new deltahat = ", fa.deltahat)
  else
    # Until we have collected at least min_updates we just sum the progress
    # values so we can later calculate their average.
    fa.deltahat += progress
  end
  fa.numupdates += 1
  if fa.numupdates == fa.min_updates
    #print("Setting deltahat, was = ", fa.deltahat, " (min_updates = $(fa.min_updates))")
    fa.deltahat = fa.deltahat / fa.min_updates
    #println(", now = ", fa.deltahat)
  end
end


text:
#! /usr/bin/julia

# Rosetta Code, Amicable pairs

function pcontrib(p::Int64, a::Int64)
    n = one(p)
    pcon = one(p)
    for i in 1:a
        n *= p
        pcon += n
    end
    return pcon
end

function divisorsum(n::Int64)
    dsum = one(n)
    for (p, a) in factor(n)
        dsum *= pcontrib(p, a)
    end
    dsum -= n
end

const L = 2*10^4
acnt = 0

println("Amicable pairs not greater than ", L)

for i in 2:L
    !isprime(i) || continue
    j = divisorsum(i)
    j < i && divisorsum(j) == i || continue
    acnt += 1
    println(@sprintf("%4d", acnt), " => ", j, ", ", i)
end

ext: jl

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
@testset "mul" begin
    rng = MersenneTwister(123456)
    P = 50
    Q = 60
    α = randn(rng)
    β = randn(rng)

    A_Matrix = randn(rng, P, Q)
    At_Matrix = collect(A_Matrix')

    A_blk_diag = BlockDiagonal([randn(rng, P, P), randn(rng, P + 1, P + 1)])
    At_blk_diag = BlockDiagonal(map(collect ∘ transpose, blocks(A_blk_diag)))

    settings = [
        (
            name="Matrix{Float64}",
            A=A_Matrix,
            At=At_Matrix,
            B=randn(rng, size(A_Matrix, 2), Q),
            C=randn(rng, size(A_Matrix, 1), Q),
        ),
        (
            name="BlockDiagonal{Float64, Matrix{Float64}}",
            A=A_blk_diag,
            At=At_blk_diag,
            B=randn(rng, size(A_blk_diag, 2), Q),
            C=randn(rng, size(A_blk_diag, 1), Q),
        ),
        (
            name="BlockDiagonal{Float64, BlockDiagonal{Float64, Matrix{Float64}}}",
            A=BlockDiagonal([A_blk_diag, A_blk_diag]),
            At=BlockDiagonal([At_blk_diag, At_blk_diag]),
            B=randn(rng, 2 * size(A_blk_diag, 2), Q),
            C=randn(rng, 2 * size(A_blk_diag, 1), Q),
        ),
    ]

    @testset "$(setting.name)" for setting in settings
        A = setting.A
        At = setting.At
        B = setting.B
        Bt = collect(B')

        C = setting.C
        Ct = collect(C')

        A_dense = collect(A)
        At_dense = collect(At)

        # Matrix-Matrix product.
        @test collect(mul!(copy(C), A, B, α, β)) ≈ α * A_dense * B + β * C
        @test collect(mul!(copy(C), At', B, α, β)) ≈ α * At_dense' * B + β * C
        @test collect(mul!(copy(C), A, Bt', α, β)) ≈ α * A_dense * Bt' + β * C
        @test collect(mul!(copy(C), At', Bt', α, β)) ≈ α * At_dense' * Bt' + β * C

        @test collect(mul!(copy(Ct), Bt, At, α, β)) ≈ α * Bt * At_dense + β * Ct
        @test collect(mul!(copy(Ct), B', At, α, β)) ≈ α * B' * At_dense + β * Ct
        @test collect(mul!(copy(Ct), Bt, A', α, β)) ≈ α * Bt * A_dense' + β * Ct
        @test collect(mul!(copy(Ct), B', A', α, β)) ≈ α * B' * A_dense' + β * Ct

        # Matrix-Vector product.
        b = B[:, 1]
        c = C[:, 1]
        @test collect(mul!(copy(c), A, b, α, β)) ≈ α * A_dense * b + β * c
    end
end


text:
using ParallelTrainingNN
using Base.Test
using Calculus
using Distributions
using Optim

@testset "Static Models" begin
    @testset "Test Linear" begin
        include("static_models/test_linear.jl")
    end
    @testset "Test Linear Sparse" begin
        include("static_models/test_linear_sparse.jl")
    end
    @testset "Test Elementwise Operator" begin
        include("static_models/test_elementwise_operator.jl")
    end
    @testset "Test Bias" begin
        include("static_models/test_bias.jl")
    end
    @testset "Test Series Model" begin
        include("static_models/test_series_model.jl")
    end
    @testset "Test Feedforward Network Model" begin
        include("static_models/test_feedforward_network.jl")
    end
end

@testset "Test Difference Equation" begin
    include("test_difference_equation.jl")
end

@testset "Test Identification Data" begin
    include("test_iddata.jl")
end

@testset "Test Levemberg-Marquadt" begin
    include("optimization/test_levenberg_marquardt.jl")
end

@testset "Test Predictor" begin
    include("test_predictor.jl")
end

@testset "Test Predictor Error" begin
    include("test_predictor_error.jl")
end

@testset "Test Preprocess Data" begin
    include("test_preprocess_data.jl")
end


text:
using ArgParse

function parse_commandline()
	s = ArgParseSettings()
	@add_arg_table! s begin
		"JSON"
      help = "path to a JSON file with the parameter settings"
      required = true
      arg_type = String
    "--output-params-dump"
      help = "path where genereated simulation parameters (SimParams) are dumped"
      arg_type = String
		"--output-run-dump-prefix"
      help = "prefix that will be added to a filename were full output for each run should be written. Each subsequent trajectory i will be saved under  \$prefix_\$i.jld2"
      arg_type = String
    "--output-daily"
      help = "path where daily trajectories should be saved"
      arg_type = String
    "--output-summary"
      help = "path where summary should be saved"
      arg_type = String
	end
  parse_args(ARGS, s)
end

max_stars_repo_name: Novirmor/MocosSimLauncher

max_issues_repo_name: Novirmor/MocosSimLauncher


text:
module MultiThreadedLU

using LinearAlgebra, LinearAlgebra.BLAS

include("hpl_seq.jl")
include("hpl_mt.jl")

export hpl_seq, hpl_mt

end


text:
# Utility function for classifier head of vision transformer-like models
seconddimmean(x) = dropdims(mean(x; dims = 2); dims = 2)

# utility function for making sure that all layers have a channel size divisible by 8
# used by MobileNet variants
function _round_channels(channels, divisor, min_value = divisor)
    new_channels = max(min_value, floor(Int, channels + divisor / 2) ÷ divisor * divisor)
    # Make sure that round down does not go down by more than 10%
    return (new_channels < 0.9 * channels) ? new_channels + divisor : new_channels
end

"""
    addrelu(x, y)

Convenience function for `(x, y) -> @. relu(x + y)`.
Useful as the `connection` argument for [`resnet`](#).
See also [`reluadd`](#).
"""
addrelu(x, y) = @. relu(x + y)

"""
    reluadd(x, y)

Convenience function for `(x, y) -> @. relu(x) + relu(y)`.
Useful as the `connection` argument for [`resnet`](#).
See also [`addrelu`](#).
"""
reluadd(x, y) = @. relu(x) + relu(y)

"""
    cat_channels(x, y, zs...)

Concatenate `x` and `y` (and any `z`s) along the channel dimension (third dimension).
Equivalent to `cat(x, y, zs...; dims=3)`.
Convenient reduction operator for use with `Parallel`.
"""
cat_channels(xy...) = cat(xy...; dims = 3)

"""
    swapdims(perm)

Convenience function for permuting the dimensions of an array.
`perm` is a vector or tuple specifying a permutation of the input dimensions.
Equivalent to `permutedims(x, perm)`.
"""
swapdims(perm) = Base.Fix2(permutedims, perm)

# Utility function for pretty printing large models
function _maybe_big_show(io, model)
    if isdefined(Flux, :_big_show)
        if isnothing(get(io, :typeinfo, nothing)) # e.g. top level in REPL
            Flux._big_show(io, model)
        else
            show(io, model)
        end
    end
end


text:
using POMDPs

#POMDPs.add("FIB")


text:
function is_leap_year(year::Int)
    rem(year,4)    == 0 &&
    (rem(year,100) != 0 ||
     rem(year,400) == 0)
end


text:
const layerincrements = (0.0, 0.025, 0.05, 0.1, 0.2, 0.3, 0.5, 1.0) .* m
const layernames = ("0", "2.5", "5", "10", "20", "30", "50", "100")
const nextlayer = 2.0m
const layerrange = (0.01, 1.20) .* m

" Load multiple years for an index "
load_point(basepath, years, shade, i::CartesianIndex, skip=()) =
    MicroclimPoint{shade, layerincrements,nextlayer,layerrange}(load_microclim(basepath, years, shade, skip, i)...)

" Load grids for a single year "
load_grid(basepath, years, shade, skip=()) =
    MicroclimGrid{shade,layerincrements,nextlayer,layerrange}(load_microclim(basepath, years, shade, skip)...)

ifhasdir(loader, basepath, file, empty) =
    if isdir(joinpath(basepath, file))
        loader()
    else
        println("$file does not exits. Skipping and returning $empty.")
        empty
    end

"""
Generic microcliate loader. An Int year and no args will return a MicroclimGrid
with full raster layers for asingle year. A range of years within 1990:2017 and a
CartesianIndex in args will return a MicroclimPoint.
"""
load_microclim(basepath, years, shade, skip, args...) = begin
    check_year(years)

    radiation = :radiation in skip ? [] : ifhasdir(basepath, "SOLR", []) do
        load_file(basepath, "SOLR", years, args...)
    end

    snowdepth = :snowdepth in skip ? [] : ifhasdir(basepath, "SNOWDEP_$(shade)pctShade", []) do
        load_file(basepath, (n, y) -> "$(n)_$(shade)pctShade/$(n)_$(shade)pctShade_$(y).nc", "SNOWDEP", years, args...)
    end

    airtemperature = :airtemperature in skip ? () : ifhasdir(basepath, "Tair120cm", ()) do
        load_file(basepath, (n, y) -> "Tair1cm_$(shade)pctShade/$(n)_$(shade)pctShade_$(y).nc", "TA1cm", years, args...),
        load_file(basepath, (n, y) -> "Tair120cm/$(n)_$(y).nc", "TA120cm", years, args...)
    end

    relhumidity = :relhumidity in skip ? () : ifhasdir(basepath, "RH120cm", ()) do
        load_file(basepath, (n, y) -> "$(n)_$(shade)pctShade/$(n)_$(shade)pctShade_$(y).nc", "RH1cm", years, args...),
        load_file(basepath, "RH120cm", years, args...)
    end

    windspeed = :windspeed in skip ? () : ifhasdir(basepath, "V120cm", ()) do
        load_file(basepath, "V1cm", years, args...), load_file(basepath, "V120cm", years, args...)
    end

    soiltemperature = :soiltemperature in skip ? () : ifhasdir(basepath, "soil$(shade)cm_$(shade)pctShade", ()) do
        load_folder(basepath, "soil", layernames, years, shade, args...)
    end

    soilwaterpotential = :soilwaterpotential in skip ? () : ifhasdir(basepath, "pot$(shade)cm_$(shade)pctShade", ()) do
        load_folder(basepath, "pot", layernames, years, shade, args...)
    end

    soilwatercontent = :soilwatercontent in skip ? () : ifhasdir(basepath, "moist$(shade)cm_$(shade)pctShade", ()) do
        load_folder(basepath, "moist", layernames, years, shade, args...)
    end

    radiation, snowdepth, airtemperature, relhumidity, windspeed,
    soiltemperature, soilwaterpotential, soilwatercontent
end

check_year(year::Integer) = year >= 1990 && year <= 2017 || error("$year is not between 1990 and 2017")
check_year(years) = check_year(years.start) && check_year(years.stop)

load_file(basepath::String, formatter, name::String, years, args...) =
    load_variable(basepath, formatter, name, name, years, args...)
load_file(basepath::String, name::String, years, args...) =
    load_variable(basepath, (n, y) -> "$(n)/$(n)_$(y).nc", name, name, years, args...)

load_folder(basepath, name, layers, years, shade, args...) =
    tuple((load_variable(basepath, (n, y) -> "$(n)$(l)cm_$(shade)pctShade/$(n)$(l)cm_$(shade)pctShade_$(y).nc", name,
                         "$(name)$(l)cm", years, args...) for l in layers)...)

"load timeseries for a single index"
load_variable(basepath::String, formatter, name::String, varname::String, years,  i::CartesianIndex) = begin
    println("Adding $varname for $years at $i")
    combined_data = Float64[]
    for year in years
        path = joinpath(basepath, formatter(name, year))
        year_data = [get_yeardata(path, varname)[i, :]...]
        combined_data = vcat(combined_data, year_data)
    end
    combined_data
end

"load timeseries for the whole grid"
load_variable(basepath::String, formatter, name::String, varname::String, years) = begin
    println("Adding $varname for $years for complete grid")

    # Initialise for the first year, as the eltype is variable
    path1 = joinpath(basepath, formatter(name, years[1]))
    years_vector = [get_yeardata(path1, varname)]

    # Add the rest of the years, if there are any
    for i in 2:length(years)
        path = joinpath(basepath, formatter(name, years[i]))
        push!(years_vector, get_yeardata(path, varname))
    end

    # Return a vector of NCDatasets for each year
    years_vector
end

function get_yeardata(path, varname)
    println("Retrieving ", varname, " from ", path)
    Dataset(path)[varname][:,:,:]
end


text:
using BufferedArrays
using Base.Test

bv = BVector{Int}(10)
@test size(bv) == (10,)
b = [i for i = 1:10]
bv .= b
@test bv == b


text:
module CrypticCrosswords

using ChartParsers
using Base.Iterators: product, drop
using Combinatorics: permutations
using DataDeps: @datadep_str
using ProgressMeter: @showprogress
import JSON

export solve, Context, IsWord, IsPrefix, IsSuffix, derive!, explain

function normalize(word)
    word |> lowercase |> strip |> (w -> replace(w, '_' => ' ')) |> (w -> replace(w, r"[^a-z0-9 ]" => ""))
end

include("SemanticSimilarity/SemanticSimilarity.jl")
using .SemanticSimilarity

include("grammar.jl")
include("solver.jl")
include("explanations.jl")
include("ptrie.jl")
include("cache.jl")

include(joinpath(@__DIR__, "..", "deps", "data_registration.jl"))
const CACHE = Cache()

is_word(x::AbstractString) = x in CACHE.words

end # module


text:
function update(s::Int,λ₀::Float64,λ₁::Float64)
	τ = rand()
	if s == 0
		if τ < λ₀
			s = 1
		end
	elseif s == 1
		if τ < λ₁
			s = 0
		end
	else
		error("What state are you in? Answer: $(s)")
	end
	return s
end


text:
module Divinity

function dat2csv(dat_path::AbstractString, csv_path::AbstractString)
    open(csv_path, write=true) do io
        for line in eachline(dat_path)
            join(io, split(line), ',')
            println(io)
        end
    end
    return csv_path
end

function dat2csv(dat_path::AbstractString)
    base, ext = splitext(dat_path)
    ext == ".dat" ||
        throw(ArgumentError("file name doesn't end with `.dat`"))
    return dat2csv(dat_path, "$base.csv")
end

end # module


text:
using Distributions
using SpecialFunctions
using JLD2
using FileIO
using Distributed

function minalpha(β)
    α =   (exp(1)* sqrt( π))/β
    return α

end

function takememory(β,v,α)
    prob =Probability_of_memory_usage(β,v,α)
    wahrscheinlikeit = rand()
    result=false
    if wahrscheinlikeit < prob
        result=true
    end
    result
end

function Probability_of_memory_usage(β,v,α)
    prob = 1- exp(-β *energymemory(v,α))
    prob
end

function energymemory(v,α)
    energy =α*exp(-(α*v)^2)*sqrt(1/ π)
    energy
end

function freeenergy(β,v,α)
    energy = log( exp(β *energymemory(v,α)) -1)/β
    energy
end

function dissipation(β,v,α)
    diss = α*α/β * (v)^2
    diss
end

function newv_drift(δ,n)
    vnew=sqrt(n)δ+randn()*δ*0.05*sqrt(n)
    vnew
end

function onestep(β,δ,oldv,α,n)

    if (takememory(β,oldv,α))
        free=freeenergy(β,oldv,α)
        diss=dissipation(β,oldv,α)
        unet=free -diss
        nextv=newv_drift(δ,n)
        n+=1
        naiveuse=0
    else
        unet = 0
        free=0
        diss=0
        nextv=newv_drift(δ,1)
        naiveuse=1
        n=2
    end

    return (unet,nextv,naiveuse,free,diss,n)
end

function livesycle(β,δ,α,meanencounters,reps)
    meanunet=0 #initializing for all the rounds for fist cost
    meanfree=0
    meandiss=0
    naiveuses = 0
    for round in 1:reps
        localnetu= 0
        netu = 0
        localfree= 0
        free = 0
        localdiss= 0
        diss = 0

        v=newv_drift(δ,1) #start the cirus at the memory wich is allways at 0
        naiveuses += 1/reps
        n=2
        for encounter in 2:meanencounters  #start at 1 because initialization is allways at 0 in the first round without any prior memory
             naiveuse =0    # Used naive for the first virus
            (netu,v,naiveuse,free,diss,n)=onestep(β,δ,v,α,n)
            localnetu+=netu
            localfree+=free
            localdiss+=diss
            naiveuses += naiveuse/reps
        end

        meanunet+= localnetu/reps
        meanfree+= localfree/reps
        meandiss+= localdiss/reps
    end

    return (meanunet,naiveuses,meanfree,meandiss)
end

function findoptalphaoflive_drift_real(β,δ,meanencounters,reps,αsteps,αmax)

    #if (minalpha(β) >αmax)
    #    println("amax is to small for the chosen beta")
    #    return 0,0,0,0,0,0
    #end
    start= 0.01
    α = range(start,stop=αmax,length=αsteps)

    data=livesycle.(β,δ,α,meanencounters,reps)
    usesmaive=zeros(αsteps)
    memuses=zeros(αsteps)
    utility=zeros(αsteps)
    free=zeros(αsteps)
    diss=zeros(αsteps)
    for ii in 1:αsteps
        utility[ii] = data[ii][1]
        usesmaive[ii] = data[ii][2]
        free[ii] = data[ii][3]
        diss[ii] = data[ii][4]
    end
    memuses= meanencounters./usesmaive

    amax= argmax(utility)
    optalpha = α[amax]
    Nusesatmax=memuses[amax]
    argmaxn = argmax(memuses)
    maxNpossible=memuses[argmaxn]
    utiliymax = utility[amax]
    freemax=free[amax]
    dissmax=diss[amax]

    return Nusesatmax,utiliymax,optalpha,maxNpossible,freemax,dissmax
end

function find_optimal_alpha_for_range_of_delta(β,meanencounters,reps,αsteps,αmax,deltasteps)

    δ = range(0.0001,stop=1.5/αmax,length= deltasteps)
    #δ = range(0.325,stop=.4,length= deltasteps)./sqrt(αmax)
    memstepsatmax =zeros(deltasteps)
    utilityopt =zeros(deltasteps)
    alphaopt =zeros(deltasteps)
    posnmax =zeros(deltasteps)
    freemax =zeros(deltasteps)
    dissmax =zeros(deltasteps)
    for i =1:deltasteps
        memstepsatmax[i],utilityopt[i],alphaopt[i],posnmax[i],freemax[i],dissmax[i]=findoptalphaoflive_drift_real(β,δ[i],meanencounters,reps,αsteps,αmax)
    end
    return memstepsatmax,utilityopt,alphaopt,posnmax,δ,freemax,dissmax

end

function find_optimal_alpha_for_range_delta_and_beta(meanencounters,reps,αsteps,αmax,deltasteps,betasteps)

    β = range(0.1,length=betasteps,stop=10.0)

    memstepsatmax= zeros(betasteps,deltasteps)
    utilityopt= zeros(betasteps,deltasteps)
    alphaopt= zeros(betasteps,deltasteps)
    posnmax= zeros(betasteps,deltasteps)
    freemax =zeros(betasteps,deltasteps)
    dissmax =zeros(betasteps,deltasteps)
    δ=zeros(deltasteps)
    Threads.@threads for ii in 1:betasteps      # for parallel execution
    # for ii in 1:betasteps                     # non parallel
    memstepsatmax[ii,:],utilityopt[ii,:],alphaopt[ii,:],posnmax[ii,:],δ,freemax[ii,:],dissmax[ii,:] = find_optimal_alpha_for_range_of_delta(β[ii],meanencounters,reps,αsteps,αmax,deltasteps)
    end
    introduction="This gives the data meory_used_at_the_max  utility_atmax  otimal_alpha posible_maximum_memoryuses delta and beta and free_at_the_max and diss_at_the_max the data is Gamma first and delta second index"
 return memstepsatmax,utilityopt,alphaopt,posnmax,δ,β,freemax,dissmax,introduction

end

max_stars_repo_name: StatPhysBio/ImmuneMemoryDM

max_issues_repo_name: StatPhysBio/ImmuneMemoryDM

max_forks_repo_name: StatPhysBio/ImmuneMemoryDM


text:
using SymbolicUtils
using StaticArrays
import Base: eltype, length, ndims, size, axes, eachindex

### Store Shape as a metadata in Term{<:AbstractArray} objects
struct ArrayShapeCtx end

#=
 There are 2 types of array terms:
   `ArrayOp{T<:AbstractArray}` and `Term{<:AbstractArray}`

 - ArrayOp represents a Einstein-notation-inspired array operation.
   it contains a field `term` which is a `Term` that represents the
   operation that resulted in the `ArrayOp`.
   I.e. will be `A*b` for the operation `(i,) => A[i,j] * b[j]` for example.
   It can be `nothing` if not known.
   - calling `shape` on an `ArrayOp` will return the shape of the array or `Unknown()`
   - do not rely on the `symtype` or `shape` information of the `.term` when looking at an `ArrayOp`.
     call `shape`, `symtype` and `ndims` directly on the `ArrayOp`.
 - `Term{<:AbstractArray}`
   - calling `shape` on it will return the shape of the array or `Unknown()`, and uses
     the `ArrayShapeCtx` metadata context to store this.

 The Array type parameter must contain the dimension.
=#

#### ArrayOp ####

"""
    ArrayOp(output_idx, expr, reduce)

A tensor expression where `output_idx` are the output indices
`expr`, is the tensor expression and `reduce` is the function
used to reduce over contracted dimensions.
"""
struct ArrayOp{T<:AbstractArray} <: Symbolic{T}
    output_idx # output indices
    expr       # Used in pattern matching
               # Useful to infer eltype
    reduce
    term
    shape
    ranges::Dict{Sym, AbstractRange} # index range each index symbol can take,
                                     # optional for each symbol
    metadata
end

function ArrayOp(T::Type, output_idx, expr, reduce, term, ranges=Dict(); metadata=nothing)
    sh = make_shape(output_idx, expr, ranges)
    ArrayOp{T}(output_idx, expr, reduce, term, sh, ranges, metadata)
end

function ArrayOp(a::AbstractArray)
    i = makesubscripts(ndims(a))
    # TODO: formalize symtype(::Type) then!
    ArrayOp(symtype(a), (i...,), a[i...], +, a)
end

ConstructionBase.constructorof(s::Type{<:ArrayOp{T}}) where {T} = ArrayOp{T}

shape(aop::ArrayOp) = aop.shape

const show_arrayop = Ref{Bool}(false)
function Base.show(io::IO, aop::ArrayOp)
    if istree(aop.term) && !show_arrayop[]
        show(io, aop.term)
    else
        print(io, "@arrayop")
        print(io, "(_[$(join(string.(aop.output_idx), ","))] := $(aop.expr))")
    end
end

symtype(a::ArrayOp{T}) where {T} = T
istree(a::ArrayOp) = true
operation(a::ArrayOp) = typeof(a)
arguments(a::ArrayOp) = [a.output_idx, a.expr, a.reduce, a.term, a.shape, a.ranges, a.metadata]

function Base.isequal(a::ArrayOp, b::ArrayOp)
    a === b && return true
    isequal(operation(a), operation(b)) &&
    isequal(a.output_idx, b.output_idx) &&
    isequal(a.expr, b.expr) &&
    isequal(a.reduce, b.reduce) &&
    isequal(a.shape, b.shape)
end

function Base.hash(a::ArrayOp, u::UInt)
    hash(a.shape, hash(a.expr, hash(a.expr, hash(a.output_idx, hash(operation(a), u)))))
end

macro arrayop(call, output_idx, expr, options...)
    @assert output_idx.head == :tuple
    rs = []
    reduce = +
    for o in options
        if isexpr(o, :call) && o.args[1] == :in
            push!(rs, :($(o.args[2]) => $(o.args[3])))
        elseif isexpr(o, :(=)) && o.args[1] == :reduce
            reduce = o.args[2]
        end
    end

    oidxs = filter(x->x isa Symbol, output_idx.args)
    iidxs = find_indices(expr)
    idxs  = union(oidxs, iidxs)
    fbody = call2term(deepcopy(expr))
    oftype(x,T) = :($x::$T)
    aop = gensym("aop")
    quote
        let
            @syms $(map(x->oftype(x, Int), idxs)...)

            expr = $fbody
            #TODO: proper Atype
            $ArrayOp(Array{$symtype(expr),
                           $(length(output_idx.args))},
                     $output_idx,
                     expr,
                     $reduce,
                     $(call2term(call)),
                     Dict($(rs...)))

        end
    end |> esc
end

const SymArray = Union{ArrayOp, Symbolic{<:AbstractArray}}
const SymMat = Union{ArrayOp{<:AbstractMatrix}, Symbolic{<:AbstractMatrix}}
const SymVec = Union{ArrayOp{<:AbstractVector}, Symbolic{<:AbstractVector}}

### Propagate ###
#
## Shape ##

function axis_in(a, b)
    first(a) >= first(b) && last(a) <= last(b)
end

function make_shape(output_idx, expr, ranges=Dict())
    matches = idx_to_axes(expr)
    for (sym, ms) in matches
        to_check = filter(m->!(shape(m.A) isa Unknown), ms)
        # Only check known dimensions. It may be "known symbolically"
        isempty(to_check) && continue
        restricted = false
        if haskey(ranges, sym)
            ref_axis = axes(ranges[sym], 1)
            restricted = true
        else
            ref_axis = axes(first(to_check).A, first(to_check).dim)
        end
        reference = axes(ref_axis, 1) # Reset to 1
        for i in (restricted ? 1 : 2):length(ms)
            m = ms[i]
            s=shape(m.A)
            if s !== Unknown()
                if restricted
                    if !axis_in(ref_axis, axes(m.A, m.dim))
                        throw(DimensionMismatch("expected $(ref_axis) to be within axes($(m.A), $(m.dim))"))
                    end
                elseif !isequal(axes(m.A, m.dim), reference)
                    throw(DimensionMismatch("expected axes($(m.A), $(m.dim)) = $(reference)"))
                end
            end
        end
    end

    sz = map(output_idx) do i
        if i isa Sym
            if haskey(ranges, i)
                return axes(ranges[i], 1)
            end
            mi = matches[i]
            @assert !isempty(mi)
            get(first(mi))
        elseif i isa Integer
            return Base.OneTo(1)
        end
    end
    # TODO: maybe we can remove this restriction?
    if any(x->x isa Unknown, sz)
        Unknown()
    else
        sz
    end
end

function ranges(a::ArrayOp)
    rs = Dict{Sym, Any}()
    ax = idx_to_axes(a.expr)
    for i in keys(ax)
        if haskey(a.ranges, i)
            rs[i] = a.ranges[i]
        else
            rs[i] = get(first(ax[i]))
        end
    end
    return rs
end

## Eltype ##

eltype(aop::ArrayOp) = symtype(aop.expr)

## Ndims ##
function ndims(aop::ArrayOp)
    length(aop.output_idx)
end

### Utils ###

# turn `f(x...)` into `term(f, x...)`
#
function call2term(expr, arrs=[])
    !(expr isa Expr) && return expr
    if expr.head == :call
        if expr.args[1] == :(:)
            return expr
        end
        return Expr(:call, term, map(call2term, expr.args)...)
    elseif expr.head == Symbol("'")
        return Expr(:call, term, adjoint, map(call2term, expr.args)...)
    end

    return Expr(expr.head, map(call2term, expr.args)...)
end

# Find all symbolic indices in expr
function find_indices(expr, idxs=[])
    !(expr isa Expr) && return idxs
    if expr.head == :ref
        return append!(idxs, filter(x->x isa Symbol, expr.args[2:end]))
    elseif expr.head == :call && expr.args[1] == :getindex || expr.args[1] == getindex
        return append!(idxs, filter(x->x isa Symbol, expr.args[3:end]))
    else
        foreach(x->find_indices(x, idxs), expr.args)
        return idxs
    end
end

struct AxisOf
    A
    dim
end

function Base.get(a::AxisOf)
    @oops shape(a.A)
    axes(a.A, a.dim)
end

function idx_to_axes(expr, dict=Dict{Sym, Vector}(), ranges=Dict())
    if istree(expr)
        if operation(expr) === (getindex)
            args = arguments(expr)
            for (axis, sym) in enumerate(@views args[2:end])
                !(sym isa Sym) && continue
                axesvec = Base.get!(() -> [], dict, sym)
                push!(axesvec, AxisOf(first(args), axis))
            end
        else
            idx_to_axes(operation(expr), dict)
            foreach(ex->idx_to_axes(ex, dict), arguments(expr))
        end
    end
    dict
end

#### Term{<:AbstractArray}
#

"""
    arrterm(f, args...; arrayop=nothing)

Create a term of `Term{<: AbstractArray}` which
is the representation of `f(args...)`.

- Calls `propagate_atype(f, args...)` to determine the
  container type, i.e. `Array` or `StaticArray` etc.
- Calls `propagate_eltype(f, args...)` to determine the
  output element type.
- Calls `propagate_ndims(f, args...)` to determine the
  output dimension.
- Calls `propagate_shape(f, args...)` to determine the
  output array shape.

`propagate_shape`, `propagate_atype`, `propagate_eltype` may
return `Unknown()` to say that the output cannot be determined

But `propagate_ndims` must work and return a non-negative integer.
"""
function arrterm(f, args...)
    atype = propagate_atype(f, args...)
    etype = propagate_eltype(f, args...)
    nd    = propagate_ndims(f, args...)

    S = if etype === Unknown() && nd === Unknown()
        atype
    elseif etype === Unknown()
        atype{T, nd} where T
    elseif nd === Unknown()
        atype{etype, N} where N
    else
        atype{etype, nd}
    end

    setmetadata(Term{S}(f, args),
                ArrayShapeCtx,
                propagate_shape(f, args...))
end

"""
    shape(s::Any)

Returns `axes(s)` or throws.
"""
shape(s) = axes(s)

"""
    shape(s::SymArray)

Extract the shape metadata from a SymArray.
If not known, returns `Unknown()`
"""
function shape(s::Symbolic{<:AbstractArray})
    if hasmetadata(s, ArrayShapeCtx)
        getmetadata(s, ArrayShapeCtx)
    else
        Unknown()
    end
end

## `propagate_` interface:
#  used in the `arrterm` construction.

atype(::Type{<:Array}) = Array
atype(::Type{<:SArray}) = SArray
atype(::Type) = AbstractArray

_propagate_atype(::Type{T}, ::Type{T}) where {T} = T
_propagate_atype(::Type{<:Array}, ::Type{<:SArray}) = Array
_propagate_atype(::Type{<:SArray}, ::Type{<:Array}) = Array
_propagate_atype(::Any, ::Any) = AbstractArray
_propagate_atype(T) = T
_propagate_atype() = AbstractArray

function propagate_atype(f, args...)
    As = [atype(symtype(T))
          for T in Iterators.filter(x->x <: Symbolic{<:AbstractArray}, typeof.(args))]
    if length(As) <= 1
        _propagate_atype(As...)
    else
        foldl(_propagate_atype, As)
    end
end

function propagate_eltype(f, args...)
    As = [eltype(symtype(T))
          for T in Iterators.filter(x->symtype(x) <: AbstractArray, args)]
    promote_type(As...)
end

function propagate_ndims(f, args...)
    error("Could not determine the output dimension of $f$args")
end

function propagate_shape(f, args...)
    error("Don't know how to propagate shape for $f$args")
end

### Wrapper type for dispatch

@symbolic_wrap struct Arr{T,N} <: AbstractArray{T, N}
    value
end

Base.hash(x::Arr, u::UInt) = hash(unwrap(x), u)

ArrayOp(x::Arr) = unwrap(x)

function Arr(x)
    A = symtype(x)
    @assert A <: AbstractArray
    Arr{maybewrap(eltype(A)), ndims(A)}(x)
end

unwrap(x::Arr) = x.value

maybewrap(T) = has_symwrapper(T) ? wrapper_type(T) : T
# These methods allow @wrapped methods to be more specific and not overwrite
# each other when defined both for matrix and vector
wrapper_type(::Type{<:AbstractMatrix}) = Arr{<:Any, 2}
wrapper_type(::Type{<:AbstractMatrix{T}}) where {T} = Arr{maybewrap(T), 2}

wrapper_type(::Type{<:AbstractVector}) = Arr{<:Any, 1}
wrapper_type(::Type{<:AbstractVector{T}}) where {T} = Arr{maybewrap(T), 1}

function Base.show(io::IO, ::MIME"text/plain", arr::Arr)
    x = unwrap(arr)
    istree(x) && print(io, "(")
    print(io, unwrap(arr))
    istree(x) && print(io, ")")
    print(io, "[", join(string.(axes(arr)), ","), "]")
end

################# Base array functions
#

# basic
# these methods are not symbolic but work if we know this info.

geteltype(s::SymArray) = geteltype(symtype(s))
geteltype(::Type{<:AbstractArray{T}}) where {T} = T
geteltype(::Type{<:AbstractArray}) = Unknown()

ndims(s::SymArray) = ndims(symtype(s))
ndims(T::Type{<:AbstractArray}) = ndims(T)

function eltype(A::Union{Arr, SymArray})
    T = geteltype(unwrap(A))
    T === Unknown() && error("eltype of $A not known")
    return T
end

function length(A::Union{Arr, SymArray})
    s = shape(unwrap(A))
    s === Unknown() && error("length of $A not known")
    return prod(length, s)
end

function size(A::Union{Arr, SymArray})
    s = shape(unwrap(A))
    s === Unknown() && error("size of $A not known")
    return length.(s)
end

function size(A::SymArray, i::Integer)
    @assert(i > 0)
    i > ndims(A) ? 1 : size(A)[i]
end

function axes(A::Union{Arr, SymArray})
    s = shape(unwrap(A))
    s === Unknown() && error("axes of $A not known")
    return s
end

function axes(A::SymArray, i)
    s = shape(A)
    s === Unknown() && error("axes of $A not known")
    return i <= length(s) ? s[i] : Base.OneTo(1)
end

function eachindex(A::Union{Arr, SymArray})
    s = shape(unwrap(A))
    s === Unknown() && error("eachindex of $A not known")
    return CartesianIndices(s)
end

function SymbolicUtils.Code.toexpr(x::ArrayOp, st)
    haskey(st.symbolify, x) && return st.symbolify[x]

    if istree(x.term)
        toexpr(x.term, st)
    else
        throw(ArgumentError("""Don't know how to turn $x
                               into code yet"""))
    end
end

function SymbolicUtils.Code.toexpr(x::Arr, st)
    toexpr(unwrap(x), st)
end

### Scalarize

scalarize(a::Array) = map(scalarize, a)
scalarize(term::Symbolic{<:AbstractArray}, idx) = term[idx...]
val2num(::Val{n}) where n = n

function replace_by_scalarizing(ex, dict)
    rule = @rule(getindex(~x, ~~i) =>
              scalarize(~x, (map(j->haskey(dict,j) ? dict[j] : j, ~~i)...,)))

    simterm = (x, f, args) -> begin
        if f === Base.literal_pow && length(args) == 3
            #=
            julia> @variables u[1:3]
            1-element Vector{Symbolics.Arr{Num, 1}}:
             u[1:3]

            julia> u.^1
            (broadcast(literal_pow, Base.RefValue{typeof(^)}(^), u, Base.RefValue{Val{1}}(Val{1}())))[1:3]
            =#
            base = args[2]
            exp = val2num(only(args[3]))
            f = only(args[1])
            f(base, exp)
        else
            f(args...)
        end
    end
    function rewrite_operation(x)
        if istree(x) && istree(operation(x))
            f = operation(x)
            replace_by_scalarizing(f, dict)(arguments(x)...)
        else
            nothing
        end
    end

    # This must be a Prewalk to avoid descending into ArrayOp (madness)
    Prewalk(Rewriters.PassThrough(Rewriters.If(x->!(x isa ArrayOp), Chain([rewrite_operation, rule]))), similarterm=simterm)(ex)
end

function scalarize(arr::AbstractArray, idx)
    arr[idx...]
end

function scalarize(arr::Term, idx)
    scalarize_op(operation(arr), arr, idx)
end

scalarize_op(f, arr) = arr

struct ScalarizeCache end

function scalarize_op(f, arr, idx)
    if hasmetadata(arr, ScalarizeCache) && getmetadata(arr, ScalarizeCache)[] !== nothing
        wrap(getmetadata(arr, ScalarizeCache)[][idx...])
    else
        thing = f(scalarize.(map(wrap, arguments(arr)))...)
        getmetadata(arr, ScalarizeCache)[] = thing
        wrap(thing[idx...])
    end
end

@wrapped function Base.:(\)(A::AbstractMatrix, b::AbstractVecOrMat)
    t = arrterm(\, A, b)
    setmetadata(t, ScalarizeCache, Ref{Any}(nothing))
end

@wrapped function Base.inv(A::AbstractMatrix)
    t = arrterm(inv, A)
    setmetadata(t, ScalarizeCache, Ref{Any}(nothing))
end

_det(x, lp) = det(x, laplace=lp)

function scalarize_op(f::typeof(_det), arr)
    det(map(wrap, collect(arguments(arr)[1])), laplace=arguments(arr)[2])
end

@wrapped function LinearAlgebra.det(x::AbstractMatrix; laplace=true)
    Term{eltype(x)}(_det, [x, laplace])
end

# A * x = b
# A ∈ R^(m x n) x ∈ R^(n, k) = b ∈ R^(m, k)
propagate_ndims(::typeof(\), A, b) = ndims(b)
propagate_ndims(::typeof(inv), A) = ndims(A)

# A(m,k) * B(k,n) = C(m,n)
# A(m,k) \ C(m,n)  = B(k,n)
function propagate_shape(::typeof(\), A, b)
    if ndims(b) == 1
        (axes(A,2),)
    else
        (axes(A,2), axes(b, 2))
    end
end

function propagate_shape(::typeof(inv), A)
    @oops shp = shape(A)
    @assert ndims(A) == 2 && reverse(shp) == shp "Inv called on a non-square matrix"
    shp
end

function scalarize(arr::ArrayOp, idx)
    @assert length(arr.output_idx) == length(idx)

    axs = ranges(arr)

    iidx = collect(keys(axs))
    contracted = setdiff(iidx, arr.output_idx)

    dict = Dict(oi => (unwrap(i) isa Symbolic ? unwrap(i) : axs[oi][i])
                for (oi, i) in zip(arr.output_idx, idx))
    partial = replace_by_scalarizing(arr.expr, dict)

    axes = [axs[c] for c in contracted]
    if isempty(contracted)
        partial
    else
        mapreduce(arr.reduce, Iterators.product(axes...)) do idx
            replace_by_scalarizing(partial, Dict(contracted .=> idx))
        end
    end
end

scalarize(arr::Arr, idx) = wrap(scalarize(unwrap(arr),
                                          unwrap.(idx)))

function scalarize(arr)
    if arr isa Arr || arr isa Symbolic{<:AbstractArray}
        map(Iterators.product(axes(arr)...)) do i
            scalarize(arr, i)
        end
    elseif istree(arr) && operation(arr) == getindex
        args = arguments(arr)
        scalarize(args[1], (args[2:end]...,))
    elseif arr isa Num
        wrap(scalarize(unwrap(arr)))
    elseif istree(arr) && symtype(arr) <: Number
        t = similarterm(arr, operation(arr), map(scalarize, arguments(arr)), symtype(arr), metadata=arr.metadata)
        scalarize_op(operation(t), t)
    else
        arr
    end
end

Base.collect(x::Arr) = scalarize(x)
Base.collect(x::SymArray) = scalarize(x)
isarraysymbolic(x) = SymbolicUtils.symtype(unwrap(x)) <: AbstractArray

Base.convert(::Type{<:Array{<:Any, N}}, arr::Arr{<:Any, N}) where {N} = scalarize(arr)


text:
#=
this is a canonical implementation
the implemented version benchmarked slightly better under a variety of inputs (julia v1.0.1)
@inline function two_sum(a::T, b::T) where {T<:AbstractFloat}
    hi = a + b
    v  = hi - a
    lo = (a - (hi - v)) + (b - v)
    return hi, lo
end
=#

const FloatWithFMA = Union{Float64, Float32, Float16}

"""
    two_sum(a, b)
Computes `hi = fl(a+b)` and `lo = err(a+b)`.
"""
@inline function two_sum(a::T, b::T) where {T<:FloatWithFMA}
    hi = a + b
    a1 = hi - b
    b1 = hi - a1
    lo = (a - a1) + (b - b1)
    return hi, lo
end

"""
    three_sum(a, b, c)

Computes `hi = fl(a+b+c)` and `md = err(a+b+c), lo = err(md)`.
"""
function three_sum(a::T, b::T, c::T) where {T<:FloatWithFMA}
    s, t   = two_sum(b, c)
    hi, u  = two_sum(a, s)
    md, lo = two_sum(u, t)
    hi, md = two_hilo_sum(hi, md)
    return hi, md, lo
end

"""
    four_sum(a, b, c, d)

Computes `hi = fl(a+b+c+d)` and `hm = err(a+b+c+d), ml = err(hm), lo = err(ml)`.
"""
function four_sum(a::T,b::T,c::T,d::T) where {T<:FloatWithFMA}
    t0, t1 = two_sum(a ,  b)
    t0, t2 = two_sum(t0,  c)
    hi, t3 = two_sum(t0,  d)
    t0, t1 = two_sum(t1, t2)
    hm, t2 = two_sum(t0, t3) # here, t0 >= t3
    ml, lo = two_sum(t1, t2)
    return hi, hm, ml, lo
end

"""
    two_sumof3(a, b, c)

Computes `hi = fl(a+b+c)` and `lo = err(a+b+c)`.
"""
function two_sumof3(a::T, b::T, c::T) where {T<:FloatWithFMA}
    s, t = two_sum(b, c)
    x, u = two_sum(a, s)
    y    = u + t
    x, y = two_sum(x, y)
    return x, y
end

@inline two_sum(a::T, b::T, c::T) where {T<:FloatWithFMA} = two_sumof3(a, b, c)

"""
    two_sumof4(a, b, c, d)

Computes `hi = fl(a+b+c+d)` and `lo = err(a+b+c+d)`.
"""
function two_sumof4(a::T, b::T, c::T, d::T) where {T<:FloatWithFMA}
    t0, t1 = two_sum(a ,  b)
    t0, t2 = two_sum(t0,  c)
    a,  t3 = two_sum(t0,  d)
    t0  = t1 + t2
    b   = t0 + t3
    return a, b
end

@inline two_sum(a::T, b::T, c::T, d::T) where {T<:FloatWithFMA} = two_sumof4(a, b, c, d)

function three_sumof4(a::T,b::T,c::T,d::T) where {T<:FloatWithFMA}
    t0, t1 = two_sum(a ,  b)
    t0, t2 = two_sum(t0,  c)
    hi, t3 = two_sum(t0,  d)
    t0, t1 = two_sum(t1, t2)
    hm, t2 = two_sum(t0, t3) # here, t0 >= t3
    ml     = t1 + t2
    return hi, hm, ml
end

@inline three_sum(a::T, b::T, c::T, d::T) where {T<:FloatWithFMA} = thre_sumof4(a, b, c, d)

"""
    two_diff(a, b)

Computes `hi = fl(a-b)` and `lo = err(a-b)`.
"""
@inline function two_diff(a::T, b::T) where {T<:FloatWithFMA}
    hi = a - b
    a1 = hi + b
    b1 = hi - a1
    lo = (a - a1) - (b + b1)
    return hi, lo
end

"""
    two_diffof3(a, b, c)

Computes `hi = fl(a-b-c)` and `lo = err(a-b-c)`.
"""
function two_diffof3(a::T, b::T, c::T) where {T<:FloatWithFMA}
    s, t = two_diff(-b, c)
    x, u = two_sum(a, s)
       y = u + t
    x, y = two_sum(x, y)
    return x, y
end

@inline two_diff(a::T, b::T, c::T) where {T<:FloatWithFMA} = two_diffof3(a, b, c)

"""
    two_hilo_sum(a, b)

*unchecked* requirement `|a| ≥ |b|`
Computes `s = fl(a+b)` and `e = err(a+b)`.
"""
@inline function two_hilo_sum(a::T, b::T) where {T<:FloatWithFMA}
    s = a + b
    e = b - (s - a)
    return s, e
end

"""
    three_hilo_sum(a, b, c)

*unchecked* requirement `|a| ≥ |b| ≥ |c|`
Computes `s = fl(a+b+c)` and `e1 = err(a+b+c), e2 = err(e1)`.
"""
function three_hilo_sum(a::T,b::T,c::T) where {T<:FloatWithFMA}
    s, t = two_hilo_sum(b, c)
    x, u = two_hilo_sum(a, s)
    y, z = two_hilo_sum(u, t)
    x, y = two_hilo_sum(x, y)
    return x, y, z
end

"""
    two_hilo_sumof3(a, b, c)

*unchecked* requirement `|a| ≥ |b| ≥ |c|`
Computes `s = fl(a+b+c)` and `e1 = err(a+b+c)`.
"""
function two_hilo_sumof3(a::T,b::T,c::T) where {T<:FloatWithFMA}
    s, t = two_hilo_sum(b, c)
    x, u = two_hilo_sum(a, s)
    y    = u + t
    x, y = two_hilo_sum(x, y)
    return x, y
end

"""
    four_hilo_sum(a, b, c, d)

*unchecked* requirement `|a| ≥ |b| ≥ |c| ≥ |d|`
Computes `s = fl(a+b+c+d)` and `e1 = err(a+b+c+d), e2 = err(e1), e3 = err(e2)`.
"""
function four_hilo_sum(a::T,b::T,c::T,d::T) where {T<:FloatWithFMA}
    t0, t1 = two_hilo_sum(a ,  b)
    t0, t2 = two_hilo_sum(t0,  c)
    hi, t3 = two_hilo_sum(t0,  d)
    t0, t1 = two_hilo_sum(t1, t2)
    hm, t2 = two_hilo_sum(t0, t3) # here, t0 >= t3
    ml, lo = two_hilo_sum(t1, t2)
    return hi, hm, ml, lo
end

"""
    three_hilo_sumof4(a, b, c, d)

*unchecked* requirement `|a| ≥ |b| ≥ |c| ≥ |d|`
Computes `s = fl(a+b+c+d)` and `e1 = err(a+b+c+d), e2 = err(e1)`.
"""
function three_hilo_sumof4(a::T,b::T,c::T,d::T) where {T<:FloatWithFMA}
    t0, t1 = two_hilo_sum(a ,  b)
    t0, t2 = two_hilo_sum(t0,  c)
    hi, t3 = two_hilo_sum(t0,  d)
    t0, t1 = two_hilo_sum(t1, t2)
    md, t2 = two_hilo_sum(t0, t3) # here, t0 >= t3
    lo     = t1 + t2
    return hi, md, lo
end

"""
    two_hilo_sumof4(a, b, c, d)

*unchecked* requirement `|a| ≥ |b| ≥ |c| ≥ |d|`
Computes `s = fl(a+b+c+d)` and `e1 = err(a+b+c+d)`.
"""
function two_hilo_sumof4(a::T,b::T,c::T,d::T) where {T<:FloatWithFMA}
    t0, t1 = two_hilo_sum(a ,  b)
    t0, t2 = two_hilo_sum(t0,  c)
    hi, t3 = two_hilo_sum(t0,  d)
    t0     = t1 + t2
    lo     = t0 + t3
    return hi, lo
end

@inline function two_prod(a::T, b::T) where {T<:FloatWithFMA}
    s = a * b
    t = fma(a, b, -s)
    return s, t
end


text:
using ChebyExp, MatrixDepot

# random 10x10
n = 1000
R = Matrix{Float64}(n,2)
A = Matrix{Float64}(10,10)
for i in 1:n
    rand!(A)
    R[i,1] = norm(expm(A)*expm(-A)-I,2)/norm(expm(A),2)
    R[i,2] = norm(chebyexp(A)*chebyexp(-A)-I,2)/norm(chebyexp(A),2)
end
a,b = mean(R,1), std(R,1)
@printf("""
Random 10x10 matrices:
expm error:     %1.2G ± %1.2G
chebyexp error: %1.2G ± %1.2G
""",a[1],b[1],a[2],b[2])

# Chebyshev and Forsythe
R = Matrix{Float64}(2,2)
for (i,A) in enumerate([matrixdepot("chebspec",10), matrixdepot("forsythe",10)])
    R[i,1] = norm(expm(A)*expm(-A)-I,2)/norm(expm(A),2)
    R[i,2] = norm(chebyexp(A)*chebyexp(-A)-I,2)/norm(chebyexp(A),2)
end
a,b = mean(R,1), std(R,1)
@printf("""
Chebyshev and Forsythe 10x10 matrices:
expm error:     %1.2G ± %1.2G
chebyexp error: %1.2G ± %1.2G
""",a[1],b[1],a[2],b[2])

# Ill-conditioned
R = Matrix{Float64}(length(matrixdepot("ill-cond")),2)
for (i,nameA) in enumerate(matrixdepot("ill-cond"))
    A = matrixdepot(nameA,Float64,16)
    norm(full(A))>5 && (A /= norm(full(A))/4)
    R[i,1] = norm(expm(full(A))*expm(-full(A))-I,2)/norm(expm(full(A)),2)
    R[i,2] = norm(chebyexp(A)*chebyexp(-A)-I,2)/norm(chebyexp(A),2)
end
a,b = mean(R,1), std(R,1)
@printf("""
Ill-conditioned 16x16 matrices:
expm error:     %1.2G ± %1.2G
chebyexp error: %1.2G ± %1.2G
""",a[1],b[1],a[2],b[2])

max_stars_repo_name: jebej/ChebyExp

max_issues_repo_name: jebej/ChebyExp

max_forks_repo_name: jebej/ChebyExp


text:
@testset "2d derivative" begin
  l = collect(-1.5:0.2:1.5);
  ins = transpose(hcat(collect.(collect(product(l, l))[:])...));
  func = (x, y) -> exp(-(x^2 + y^2));
  dfunc = (x, y) -> 4 * x * y * exp(-(x^2 + y^2));
  fs = [func(ins[i, 1], ins[i, 2]) for i in 1:size(ins, 1)];
  dfs = [dfunc(ins[i, 1], ins[i, 2]) for i in 1:size(ins, 1)];

  m = mwls_kd(hcat(ins, fs), 1, (d, eps) -> exp(-d^2), maxDegree = 5)
  cll = mwls_cll(ins, fs, 1, (d, eps) -> exp(-d^2), maxDegree = 5)

  teststatus = true
  for i in 1:size(ins, 1)
    teststatus &= isapprox(mwls_diff(m, [ins[i, 1], ins[i, 2]], (1, 1)),
                           dfunc(ins[i, 1], ins[i, 2]),
                           atol = 0.1)
    if !teststatus
      @show ins[i, :]
      @show mwls_diff(m, [ins[i, 1], ins[i, 2]], (1, 1))
      @show dfunc(ins[i, 1], ins[i, 2])
      break
    end
  end
  @test teststatus == true

  teststatus = true
  for i in 1:size(ins, 1)
    teststatus &= isapprox(mwls_diff(cll, [ins[i, 1], ins[i, 2]], (1, 1)),
                           dfunc(ins[i, 1], ins[i, 2]),
                           atol = 0.1)
    if !teststatus
      @show ins[i, :]
      @show mwls_diff(m, [ins[i, 1], ins[i, 2]], (1, 1))
      @show dfunc(ins[i, 1], ins[i, 2])
      break
    end
  end
  @test teststatus == true
end


text:
// Variable shadow in method defintion.

int main () {
  Counter c;
  c = new Counter;
  c.incr(1);
  c.incr(20);
  c.incr(100);
  int x = c.value();
  printInt(x);
  return 0;
}

class Counter {
  int val;

  void incr (int val) {val = val; return;}

  int value () {return val;}

}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
function derivative(f, x::Union{Number,AbstractArray{<:Number}},
                    integrator)
    local d
    alg = unwrap_alg(integrator, true)
    if get_current_alg_autodiff(integrator.alg, integrator.cache)
      d = ForwardDiff.derivative(f, x)
    else
      d = DiffEqDiffTools.finite_difference_gradient(f, x, alg.diff_type, eltype(x), Val(false))
    end
    d
end

function derivative!(df::AbstractArray{<:Number}, f, x::Union{Number,AbstractArray{<:Number}}, fx::AbstractArray{<:Number}, integrator::DEIntegrator)
    if alg_autodiff(integrator.alg)
        ForwardDiff.derivative!(df, f, fx, x)
    else
        RealOrComplex = eltype(integrator.u) <: Complex ? Val{:Complex} : Val{:Real}
        DiffEqDiffTools.finite_difference!(df, f, x, integrator.alg.diff_type, RealOrComplex, fx)
    end
    nothing
end

jacobian_autodiff(f,x,_,_)=ForwardDiff.derivative(f,x)
jacobian_autodiff(f,x::AbstractArray,_,hascolorvec::Val{false})=ForwardDiff.jacobian(f,x)
function jacobian_autodiff(f,x::AbstractArray,integrator,hascolorvec::Val{true})
  colorvec=integrator.f.colorvec
  jac=integrator.f.jac_prototype
  J=jac isa SparseMatrixCSC ? similar(jac) : fill(0.,size(jac))
  forwarddiff_color_jacobian!(J,f,x,colorvec=colorvec,sparsity=jac)
  J
end
jacobian_finitediff(f,x,difftype,_,_)=DiffEqDiffTools.finite_difference_derivative(f, x, difftype, eltype(x))
jacobian_finitediff(f,x::AbstractArray,difftype,_,hascolorvec::Val{false})=DiffEqDiffTools.finite_difference_jacobian(f, x, difftype, eltype(x))
jacobian_finitediff(f,x::AbstractArray,difftype,integrator,hascolorvec::Val{true})=
  DiffEqDiffTools.finite_difference_jacobian(f, x, difftype, eltype(x), colorvec=integrator.f.colorvec,sparsity=integrator.f.jac_prototype)

function jacobian(f, x,
                  integrator::DiffEqBase.DEIntegrator)
    local J
    alg = unwrap_alg(integrator, true)
    if get_current_alg_autodiff(integrator.alg, integrator.cache)
      J = jacobian_autodiff(f,x,integrator,Val(DiffEqBase.has_colorvec(integrator.f)))
    else
      J = jacobian_finitediff(f,x,alg.diff_type,integrator,Val(DiffEqBase.has_colorvec(integrator.f)))
    end
    J
end

jacobian_autodiff!(J,f,fx,x,jac_config::ForwardColorJacCache)=forwarddiff_color_jacobian!(J,f,x,jac_config)
jacobian_autodiff!(J,f,fx,x,jac_config::ForwardDiff.JacobianConfig)=ForwardDiff.jacobian!(J, f, fx, x, jac_config)

function jacobian!(J::AbstractMatrix{<:Number}, f, x::AbstractArray{<:Number}, fx::AbstractArray{<:Number}, integrator::DEIntegrator, jac_config)
    if alg_autodiff(integrator.alg)
      jacobian_autodiff!(J, f, fx, x, jac_config)
    else
      DiffEqDiffTools.finite_difference_jacobian!(J, f, x, jac_config)
    end
    nothing
end

jac_cache_autodiff(alg,f,uf,du1,uprev,u,hascolorvec::Val{false})=ForwardDiff.JacobianConfig(uf,du1,uprev,ForwardDiff.Chunk{determine_chunksize(u,alg)}())
jac_cache_autodiff(alg,f,uf,du1,uprev,u,hascolorvec::Val{true})=ForwardColorJacCache(uf,uprev,colorvec=f.colorvec,sparsity=f.jac_prototype)

function DiffEqBase.build_jac_config(alg::StochasticDiffEqAlgorithm,f,uf,du1,uprev,u,tmp,du2)
  if !has_jac(f)
    if alg_autodiff(alg)
      jac_config = jac_cache_autodiff(alg,f,uf,du1,uprev,u,Val(DiffEqBase.has_colorvec(f)))
    else
      colorvec= f.colorvec isa Nothing ? Base.OneTo(length(u)) : f.colorvec
      if alg.diff_type != Val{:complex}
        jac_config = DiffEqDiffTools.JacobianCache(tmp,du1,du2,alg.diff_type,colorvec=colorvec,sparsity=f.jac_prototype)
      else
        jac_config = DiffEqDiffTools.JacobianCache(Complex{eltype(tmp)}.(tmp),Complex{eltype(du1)}.(du1),nothing,alg.diff_type,colorvec=colorvec,sparsity=f.jac_prototype)
      end
    end
  else
    jac_config = nothing
  end
  jac_config
end


text:
struct ZygoteRuleConfig{CTX<:AContext} <: RuleConfig{Union{HasReverseMode,NoForwardsMode}}
  context::CTX
end
ZygoteRuleConfig() = ZygoteRuleConfig(Context())

const rrule_fallback_method = Base.which(rrule, Tuple{Any, Vararg{Any}})
const rrule_redispatcher_method = Base.which(rrule, Tuple{RuleConfig, Any, Vararg{Any}})

"""
  has_chain_rrule(T)

For a type-tuple `T` e.g. `Tuple{typeof(f), Int, Float64}`, checks if there is a `rrule` defined for it.
Excluding the generic fallback.
The first return value is `true` if the `rrule` exists, `false` otherwise.
If it does not, then the second argument is a list of edges to attach to the CodeInfo for a generated function,
such that if a suitable rule is defined later, the generated function will recompile.
"""
function has_chain_rrule(T)
  config_T, arg_Ts = Iterators.peel(T.parameters)
  m_with_config = meta(Tuple{typeof(rrule), config_T, arg_Ts...})
  if m_with_config.method === rrule_redispatcher_method
    # it is being redispatched without config, so check it that hits the fallback
    m_without_config = meta(Tuple{typeof(rrule), arg_Ts...})
    if m_without_config.method === rrule_fallback_method
      # no rrule exists, return instance for m_with_config as that will be invalidated
      # directly if configured rule added, or indirectly if unconfigured rule added
      return false, m_with_config.instance
    end
  end
  # otherwise found a rrule, no need to add any edges, as it will generate code with
  # natural edges.
  return true, nothing
end

"""
    is_kwfunc(sigt...)

Determines if `sigt` is the type signature of a kwfunction.
Each element of `sigt` should be a type.
Either the first 3 types are a kwfunc type, a NamedTuple and the matching base function type,
or the first argument is the base function type and it is not a kwfunction.
the remaining types in `sigt` are the types of the argument.

"""
is_kwfunc(::Vararg{Any}) = false
is_kwfunc(k, ::Type{<:NamedTuple}, f, args...) = k===Core.kwftype(f)

"""
    wrap_chainrules_output(x)

Convert `x` from the differentials types ChainRules uses to the format Zygote uses internally.
"""
@inline wrap_chainrules_output(x) = unthunk(x)  # For now we are just not going to deal with thunks
@inline wrap_chainrules_output(x::Tuple) = map(wrap_chainrules_output, x)
# Zygote convention: even if many AbstractZero partials (i.e. multi-input function), make just 1 nothing.
@inline wrap_chainrules_output(x::Tuple{Vararg{ChainRules.AbstractZero}}) = nothing
@inline wrap_chainrules_output(x::ChainRules.AbstractZero) = nothing
for T_outer in (:Tuple, :NamedTuple)
  # we create separate methods rather than using a `Union` + an `if` so that we avoid a
  # branch that changes output type, because nested AD on that kinda thing makes Zygote less
  # than happy.
  @eval @inline function wrap_chainrules_output(x::ChainRules.Tangent{P, T}) where {P, T<:$T_outer}
    xp = map(wrap_chainrules_output, canonicalize(x))
    convert($T_outer, xp)
  end
end

"""
    wrap_chainrules_input(x)

Convert `x` from the format Zygote uses internally to differentials types ChainRules uses.
"""
@inline wrap_chainrules_input(x) = x
@inline wrap_chainrules_input(::Nothing) = ChainRules.ZeroTangent()
@inline function wrap_chainrules_input(xs::Union{Tuple, NamedTuple})
  xp = map(wrap_chainrules_input, xs)
  ChainRules.Tangent{Any, typeof(xp)}(xp)
end

"""
  ZBack{F}(back) <: Function

Wrapper for a ChainRules pullback `back`, that causes it to follow Zygote conventions.
(A functor here is used rather than a closure to avoid boxing issues);
"""
struct ZBack{F} <: Function
  back::F
end
@inline (s::ZBack)(dy) = wrap_chainrules_output(s.back(wrap_chainrules_input(dy)))
# `nothing->nothing` can be deleted after https://github.com/FluxML/Zygote.jl/issues/603
# though it might be worth keeping as a performance optimization (benchmarking pending)
@inline (s::ZBack)(::Nothing) = nothing

"""
    chain_rrule(config, f, args...)

Returns a the (primal) value of `f(args...)` and a pullback, by invoking `ChainRulesCore.rrule(f, args...)`.
The pullback is appropriately wrapped up to follow Zygote conventions.
"""
@inline function chain_rrule(config, f, args...)
  y, back = rrule(config, f, args...)
  return y, ZBack(back)
end

"""
  chain_rrule_kw(config, kwf, kwargs, f, args...)

As per [`chain_rrule`](@ref) but with support for kwargs.
`kwf` should be the kwfunc matching to `f`, and `kwargs` are a `NamedTuple` of keyword arguments.
"""
@inline function chain_rrule_kw(config, kwf, kwargs, f, args...)
  y, back = rrule(config, f, args...; kwargs...)
  function kw_zpullback(dy)
    dxs = ZBack(back)(dy)
    if dxs === nothing  # if dxs is nothing, then all partiaols are nothing
      # Zygote convention is a single nothing no mather how partials, if all are nothing
      return nothing
    else
      return (nothing, nothing, dxs...)  # first two nothings are for kwfunc and kwargs
    end
  end
  return y, kw_zpullback
end

function ChainRulesCore.rrule_via_ad(config::ZygoteRuleConfig, f, args...)
    y, pb = _pullback(config.context, f, args...)
    ad_pullback(Δ) = zygote2differential(pb(wrap_chainrules_output(Δ)), (f, args...))
    return y, ad_pullback
end

"""
    zygote2differential(x)

Convert input `x` from the Zygote format to the ChainRules differential types.
"""
zygote2differential(x, primal) = z2d(x, primal)
zygote2differential(::Nothing, ::Any) = NoTangent()
zygote2differential(t::Tuple, primal::Tuple) = map(z2d, t, primal)
zygote2differential(t::Tuple, primal) = (@warn "primal should be a tuple, not $primal"; return t)
z2d(x, ::Any) = x
z2d(::Nothing, ::Any) = NoTangent()
z2d(a::AbstractArray{<:Number}, primal::AbstractArray{T}) where T = a
z2d(a::AbstractArray, primal::AbstractArray{T}) where T = z2d.(a, primal)
z2d(x::Union{AbstractZero, Tangent}, ::Any) = (difftype_warn(x); return x)
function z2d(t::Tuple, primal::Tuple)
  tp::Tuple = map(z2d, t, primal)
  primal_type = typeof(primal)
  return canonicalize(Tangent{primal_type, typeof(tp)}(tp))
end

function z2d(t::NamedTuple, primal)
  primal_type = typeof(primal)
  fnames = fieldnames(primal_type)
  complete_t = NamedTuple{fnames}(fn in keys(t) ? t[fn] : nothing for fn in fnames)
  primals = NamedTuple{fnames}(getfield(primal, fn) for fn in fnames)
  tp::NamedTuple = map(z2d, complete_t, primals)
  return canonicalize(Tangent{primal_type, typeof(tp)}(tp))
end


text:
"""
    Crash

Collection of crash methods used to generate initial decisions in structured algorithms.

...
# Available crash methods
- [`None`](@ref)
- [`EVP`](@ref)
- [`Scenario`](@ref)
- [`Custom`](@ref)
...

## Examples

The following solves a stochastic program `sp` created in `StochasticPrograms.jl` using an L-shaped algorithm with trust-region and Clp as an `lpsolver` and by generating an initial decision with the `EVP` crash.

```jldoctest
julia> optimize!(sp, solver = LShapedSolver(GLPKSolverLP(), crash=Crash.EVP(), regularize = TrustRegion()))
L-Shaped Gap  Time: 0:00:00 (8 iterations)
  Objective:       -855.8333333333339
  Gap:             0.0
  Number of cuts:  4
  Iterations:      8
:Optimal
```
"""
module Crash

using StochasticPrograms
using MathProgBase

abstract type CrashMethod end

"""
    None

Randomize the initial decision (default).

"""
struct None <: CrashMethod end

function (::None)(stochasticprogram::StochasticProgram, solver::MathProgBase.AbstractMathProgSolver)
    return rand(decision_length(stochasticprogram))
end

"""
    EVP

Solve the expected value problem corresponding to the stochastic program and use the expected value solution as initial decision.

"""
struct EVP <: CrashMethod end

function (::EVP)(sp::StochasticProgram, solver::MathProgBase.AbstractMathProgSolver)
    evp = StochasticPrograms.EVP(sp; solver = solver)
    status = solve(evp)
    status != :Optimal && error("Could not solve EVP model during crash procedure. Aborting.")
    return evp.colVal[1:decision_length(sp)]
end

"""
    Scenario

Solve the wait-and-see problem corresponding a supplied scenario and use the optimal solution as initial decision.

"""
struct Scenario{S <: AbstractScenario} <: CrashMethod
    scenario::S

    function Scenario(scenario::S) where S <: AbstractScenario
        return new{S}(scenario)
    end
end

function (crash::Scenario)(so::StochasticProgram, solver::MathProgBase.AbstractMathProgSolver)
    ws = WS(sp, crash.scenario; solver = solver)
    status = solve(ws)
    status != :Optimal && error("Could not solve wait-and-see model during crash procedure. Aborting.")
    return ws.colVal[1:decision_length(sp)]
end

"""
    Custom(x₀)

Use the user-supplied `x₀` as initial decision.

"""
struct Custom{T <: AbstractFloat} <: CrashMethod
    x₀::Vector{T}

    function Custom(x₀::Vector{T}) where T <: AbstractFloat
        return new{T}(x₀)
    end
end

function (crash::Custom)(sp::StochasticProgram, solver::MathProgBase.AbstractMathProgSolver)
    return crash.x₀[1:decision_length(sp)]
end

end

CrashMethod = Crash.CrashMethod


text:
# -------------------- Loading and Saving Image Utilities --------------------

im_mean = reshape([0.485, 0.456, 0.406], 1, 1, 3) |> gpu
im_std = reshape([0.229, 0.224, 0.225], 1, 1, 3) |> gpu

"""
    load_image(path, resize = false; size_save = true)

Load the original image as an array.

Arguments:
1. `path`: Path to the original image.
2. `resize`: Defaults to `false`. Setting it to `true` will reshape
             the original image to 224×224. This is recommended when
             using on `cpu`. However, when the image is finally saved
             the original dimensions are restored.
3. `size_save`: Set it to `false` when the size is not to be stored
                globally.
"""

function load_image(path, resize = false; size_save = true)
    img = load(path)
    if(size_save)
        global original_size = size(img)
    end
    if(resize)
        image_to_arr(imresize(img, (224, 224)))
    else
        image_to_arr(img)
    end
end

"""
    load_guide_image(path)

Utility to load the guide image. This essentially falls back to the
`load_image` function. The guide image is always resized to 224×224.
Simce most of the models used are trained on `Imagenet` the features are
best captured when the dimensions are 224×224.

Arguments:
1. `path`: Path to the guiding image. Ideally the guiding image should have
           a dimension close to 224×224 to ensure not a lot of features are
           lost while downscaling.
"""

load_guide_image(path) = load_image(path, true, size_save = false)

"""
    generate_image(x, resize_original = false)

Get back the original image from its array. It is assumed that the image was
preprocessed and is an RGB image.

Arguments:
1. `x`: Array from which the image must be generated.
2. `resize_original`: Set it to `true` if the original dimensions of the image
                      were changed
"""

function generate_image(x, resize_original = false)
    x = reshape(x, size(x)[1:3]...)
    x = x .* im_std .+ im_mean
    x = clamp.(permutedims(x, [3,2,1]), 0, 1) |> cpu
    if resize_original
        imresize(colorview(RGB, x), original_size)
    else
        colorview(RGB, x)
    end
end

"""
    save_image(path, x)

Function to save the generated image. Internally it simply makes a call to
the `generate_image` function.

Arguments:
1. `path`: The path where the image must be saved.
2. `x`: The array of the image to be saved.
"""

save_image(path, x) = save(path, generate_image(x, true))

"""
    image_to_arr(img; preprocess = true)

Converts an image to an array. Also applies preprocessing using the Imagenet
mean and standard deviation. The eltype of the array is always of type Float32.

Arguments:
1. `img`: Image whose array is to be generated.
2. `preprocess`: If set to `true` the image is preprocessed wrt Imagenet.
"""

function image_to_arr(img; preprocess = true)
    local x = img
    x = Float32.(channelview(img))
    x = permutedims(x, [3,2,1]) |> gpu
    if(preprocess)
        x = (x .- im_mean)./im_std
    end
    x = reshape(x, size(x,1), size(x,2), size(x,3), 1)
end

# -------------------- Utilities to modify Images --------------------

"""
    zoom_image(x, scale_x, scale_y)

Function to centre zoom into the image.

Arguments:
1. `x`: The array of the image.
2. `scale_x`: Scale along the x axis.
3. `scale_y`: Scale along the y axis.
"""

function zoom_image(x, scale_x, scale_y)
    img = generate_image(x)
    inv_scalex = 1 / scale_x
    inv_scaley = 1 / scale_y
    local s = size(img)
    local r1 = clamp(Int(ceil(s[1] * (1 - inv_scalex) / 2)), 1, s[1]) : clamp(Int(ceil(s[1] * (1 + inv_scalex) / 2)), 1, s[1])
    local r2 = clamp(Int(ceil(s[2] * (1 - inv_scaley) / 2)), 1, s[2]) : clamp(Int(ceil(s[2] * (1 + inv_scaley) / 2)), 1, s[2])
    img = imresize(img[r1, r2], s)
    image_to_arr(img)
end

# -------------------- Loading the Neural Network Model --------------------

"""
    load_model(layer, m = VGG19)

Function that globally loads the model.

Arguments:
1. `layer`: The layer upto which the model is to be imported.
2. `m`: The model to be loaded. All `Metalhead.jl` models are supported by
        default. To use a custom model make sure that the layers of the model
        are indexable. For more information about this look at Metalhead.jl.
"""

function load_model(layer, m = VGG19)
    model = m()
    global model = Chain(model.layers[1:layer]...) |> gpu
end

# ---------------- Utilities for Calculations on Matrices -------------------

"""
    argmax(A, dims)

Find the maximum index along a particular dimension

Arguments:
1. `A`: The matrix
2. `dims`: The dimension along which the maxima needs to be found.
"""

function argmax(A, dims)
   z = findmax(A, dims)[2] .% size(A, dims)
   z[z.==0] .= size(A,dims)
   z
end

#-----------------Utilities to apply Deepdream on a Batch--------------------

"""
    dream_batch(depth, iterations, η, octave_scale, num_octaves, guide = 1.0; path = "../images/")

Utility to generate dreams on a batch of images. This is a recommended method to be used
when applying deepdream on a bunch of images without varying the hyperparameters.

For individual parameter description look into the documentation of `deepdream` function.
"""

function dream_batch(depth, iterations, η, octave_scale, num_octaves, guide = 1.0; path = "../images/")
    images = [joinpath(path, i) for i in readdir(path)]
    load_model(depth)
    for i in images
        save_path = *((split(i, ".")[1:end-1])...) * "_dream.jpg"
        deepdream(load_image(i), iterations, η, octave_scale, num_octaves, save_path, guide)
        info("Image saved at $(save_path)")
    end
end

#------------------Video Generation Capabilities-----------------------------

"""
    writevideo(fname; overwrite=true, fps=1, options=``)

Must be called only after the `recurdream` function is called. This generates
a video using the frames generated. To use this `ffmpeg` must be present on
the system.

Arguments:
1. `fname`: Path of the video to be generated.

Other arguments are specific to `ffmpeg`. If a lot of frames are available, the
fps should be increased to generate better quality videos.
"""

function writevideo(fname; overwrite=true, fps=1, options=``)
    ow = overwrite ? `-y` : `-n`
    nframes = length(frames)
    h, w = size(frames[1])
    open(`ffmpeg
            -loglevel warning
            $ow
            -f rawvideo
            -pix_fmt rgb24
            -s:v $(h)x$(w)
            -r $fps
            -i pipe:0
            $options
            -vf "transpose=0"
            -pix_fmt yuv420p
            $fname`, "w") do out
        for i = 1:nframes
            write(out, frames[i])
        end
    end
end


text:
# Test I/O
module Mock
type MockGenerator <: Component
    model::EnergyModel
end

@testset "NetCDF" begin
    model = EnergyModel()
    MockGenerator(model)
    naming(model, :)
    defVar()

    io = EM.NcIO("test.nc")
    @test EM.get(io, :generator) == nothing
end

max_stars_repo_name: coroa/EnergyModels

max_issues_repo_name: coroa/EnergyModels

max_forks_repo_name: coroa/EnergyModels


text:
"""
General framework for representing Google JSON APIs.
"""
module api

export APIRoot, APIResource, APIMethod, set_session!, get_session, iserror

using Dates, Base64

using HTTP
import MbedTLS
import Libz
import JSON
using Markdown

using ..session
using ..error
using ..root

const GZIP_MAGIC_NUMBER = UInt8[0x1f, 0x8b, 0x08]

"""
    path_tokens(path)

Extract tokens from a path, e.g.

```
path_tokens("/{foo}/{bar}/x/{baz}")
# output
3-element Array{SubString{String},1}:
 "{foo}"
 "{bar}"
 "{baz}"
```
"""
path_tokens(path::AbstractString) = collect((m.match for m = eachmatch(r"{\w+}", path)))

"""
    path_replace(path, values)

Replace path tokens in path with values.

Assumes values are provided in the same order in which tokens appear in the path.

```
path_replace("/{foo}/{bar}/{baz}", ["this", "is", "it"])
# output
"/this/is/it"
```
"""
path_replace(path::AbstractString, values) = reduce(
                (x, y) -> replace(x, y[1]=>HTTP.URIs.escapeuri(y[2]), count=1),
                zip(path_tokens(path), values); init=path)
#function path_replace(path::AbstractString, values)
#    for value in values
#        path = replace(path, r"{\w+}"=>value, count=1)
#    end
#    path
#    #reduce((x, y) -> replace(x, y[1], HTTP.URIs.escapeuri(y[2]), 1), path,
#    #       zip(path_tokens(path), values))
#end

"""Check if response is/contains an error"""
iserror(x::AbstractDict{Symbol}) = haskey(x, :error)
iserror(::Any) = false

"""
    APIMethod(verb, path, description)

Maps a method in the API to an HTTP verb and path.
"""
struct APIMethod
    verb::Symbol
    path::String
    description::String
    default_params::Dict{Symbol, Any}
    transform::Function
    function APIMethod(verb::Symbol, path::AbstractString, description::AbstractString,
                       default_params::AbstractDict{Symbol}=Dict{Symbol, Any}();
                       transform=(x, t) -> x)
        new(verb, path, description, default_params, transform)
    end
end
function Base.print(io::IO, x::APIMethod)
    println(io, "$(x.verb): $(x.path)")
    Markdown.print_wrapped(io, x.description)
end
Base.show(io::IO, x::APIMethod) = print(io, x)

"""
    APIResource(path, methods)

Represents a resource in the API, typically rooted at a specific point in the
REST hierarchy.
"""
struct APIResource
    path::String
    methods::Dict{Symbol, APIMethod}
    transform::Union{Function, DataType}
    function APIResource(path::AbstractString, transform=identity; methods...)
        if isempty(methods)
            throw(APIError("Resource must have at least one method."))
        end
        methods = Dict(methods)
        # build out non-absolute paths
        if isurl(path)
            for (name, method) in methods
                if isempty(method.path)
                    method_path = path
                elseif !isurl(method.path)
                    method_path = startswith(method.path, ":") ? "$(path)$(method.path)" : "$(path)/$(method.path)"
                else
                    continue
                end
                methods[name] = APIMethod(method.verb, method_path, method.description, method.default_params; transform=method.transform)
            end
        end
        new(path, methods, transform)
    end
end
function Base.print(io::IO, x::APIResource)
    println(io, x.path, "\n")
    for (name, method) in sort(collect(x.methods), by=(x) -> x[1])
        println(io)
        println(io, name)
        println(io, repeat("-", length(string(name))))
        print(io, method)
    end
end
Base.show(io::IO, x::APIResource) = print(io, x)

"""
    APIRoot(...)

Represent a Google JSON API containing resources, accessible via scopes.
"""
struct APIRoot
    path::String
    scopes::Dict{String, String}
    resources::Dict{Symbol, APIResource}
    """
        APIRoot(path, scopes; resources...)

    An API rooted at `path` with specified OAuth 2.0 access scopes and
    resources.
    """
    function APIRoot(path::AbstractString, scopes::AbstractDict{<: AbstractString, <: AbstractString}; resources...)
        if !isurl(path)
            throw(APIError("API root must be a valid URL."))
        end
        if isempty(resources)
            throw(APIError("API must contain at least one resource."))
        end
        resources = Dict(resources)
        # build out non-absolute paths
        for (name, resource) in resources
            if isempty(resource.path)
                resource_path = path
            elseif !isurl(resource.path)
                resource_path = "$(path)/$(resource.path)"
            else
                continue
            end
            resources[name] = APIResource(resource_path, resource.transform; resource.methods...)
        end
        new(path, scopes, resources)
    end
end

function Base.print(io::IO, x::APIRoot)
    println(io, x.path, "\n")
    for (name, resource) in sort(collect(x.resources), by=(x) -> x[1])
        println(io)
        println(io, name)
        println(io, repeat("=", length(string(name))))
        println(io, resource)
        println(io, repeat("-", 79))
    end
end
Base.show(io::IO, x::APIRoot) = print(io, x)

"""
    set_session!(api, session)

Set the default session for a specific API. Set session to `nothing` to
forget session.
"""
function set_session!(api::APIRoot, session::Union{GoogleSession, Nothing})
    _default_session[api] = session
    nothing
end

"""
    get_session(api)

Get the default session for a specific API.
"""
function get_session(api::APIRoot)
    get(_default_session, api, nothing)
end

function (api::APIRoot)(resource_name::Symbol)
    resource = try api.resources[resource_name] catch
        throw(APIError("Unknown resource type: $resource_name"))
    end
end

function (api::APIRoot)(resource_name::Symbol, method_name::Symbol, args...; kwargs...)
    resource = api(resource_name)
    method = try resource.methods[method_name] catch
        throw(APIError("Unknown method for resource $resource_name: $method_name"))
    end
    kwargs = Dict(kwargs)
    session = pop!(kwargs, :session, get(_default_session, api, nothing))
    if session === nothing
        throw(SessionError("Cannot use API without a session."))
    end
    execute(session, resource, method, args...; kwargs...)
end

"""
    execute(session::GoogleSession, resource::APIResource, method::APIMethod,
            path_args::AbstractString...[; ...])

Execute a method against the provided path arguments.

Optionally provide parameters and data (with optional MIME content-type).
"""
function execute(session::GoogleSession, resource::APIResource, method::APIMethod,
            path_args::AbstractString...;
            data::Union{AbstractString, AbstractDict, Vector{UInt8}}=HTTP.nobody,
            gzip::Bool=false, content_type::AbstractString="application/json",
            debug::Bool=false, raw::Bool=false,
            max_backoff::TimePeriod=Second(64), max_attempts::Int64=10,
            params...)
    if length(path_args) != path_tokens(method.path) |> length
        throw(APIError("Number of path arguments do not match"))
    end

    # obtain and use access token
    auth = authorize(session)
    headers = Dict{String, String}(
        "Authorization" => "$(auth[:token_type]) $(auth[:access_token])"
    )
    params = Dict(params)

    # check if data provided when not expected
    if xor((!isempty(data)), in(method.verb, (:POST, :UPDATE, :PATCH, :PUT)))
        data = HTTP.nobody
        content_type = ""
        headers["Content-Length"] = "0"
    end

    if !isempty(content_type)
        headers["Content-Type"] = content_type
    end

    # serialise data to JSON if necessary
    if !isempty(data)
        if isa(data, AbstractDict) || content_type == "application/json"
            data = JSON.json(data)
        elseif isempty(data)
            headers["Content-Length"] = "0"
        end
        if gzip
            params[:contentEncoding] = "gzip"
            if !all(data[1:3] .== GZIP_MAGIC_NUMBER)
                # check the data compression using gzip magic number
                data = read(Vector{UInt8}(data) |> Libz.ZlibDeflateInputStream)
            end
        end
    end

    # merge in default parameters and evaluate any expressions
    params = merge!(copy(method.default_params), Dict(params))
    extra = Dict(:project_id => session.credentials.project_id)
    for (key, val) in params
        if isa(val, Symbol)
            params[key] = extra[val]
        end
    end

    # attempt request until exceeding maximum attempts, backing off exponentially
    res = nothing
    max_backoff = Millisecond(max_backoff)
    for attempt = 1:max(max_attempts, 1)
        if debug
            @info("Attempt: $attempt")
        end
        res = try
            HTTP.request(string(method.verb),
                        path_replace(method.path, path_args), headers, data;
                        query=params )
        catch e
        #    if isa(e, Base.IOError) &&
        #        e.code in (Base.UV_ECONNRESET, Base.UV_ECONNREFUSED, Base.UV_ECONNABORTED,
        #                   Base.UV_EPIPE, Base.UV_ETIMEDOUT)
        #    elseif isa(e, MbedTLS.MbedException) &&
        #            e.ret in (MbedTLS.MBEDTLS_ERR_SSL_TIMEOUT, MbedTLS.MBEDTLS_ERR_SSL_CONN_EOF)
        #    else
        #        println("get a HTTP request error: ", e)
        #        rethrow(e)
        #    end
            println("get a HTTP request error: ", e)
            rethrow(e)
        end

        if debug && (res !== nothing)
            @info("Request URL: $(res.request.target)")
            @info("Response Headers:\n" * join(("  $name: $value" for (name, value) in
                                                sort(collect(res.headers))), "\n"))
            @info("Response Data:\n  " * base64encode(res.body))
            @info("Status: ", res.status)
        end

        # https://cloud.google.com/storage/docs/exponential-backoff
        if (res === nothing) || (div(res.status, 100) == 5) || (res.status == 429)
            if attempt < max_attempts
                backoff = min(Millisecond(floor(Int, 1000 * (2 ^ (attempt - 1) + rand()))), max_backoff)
                warn("Unable to complete request: Retrying ($attempt/$max_attempts) in $backoff")
                sleep(backoff / Millisecond(Second(1)))
            else
                warn("Unable to complete request: Stopping ($attempt/$max_attempts)")
            end
        else
            break
        end
    end

    # if response is JSON, parse and return. otherwise, just dump data
    # HTTP response header type is Vector{Pair{String,String}}
    # https://github.com/JuliaWeb/HTTP.jl/blob/master/src/Messages.jl#L166
    for (key, value) in res.headers
        if key=="Content-Type"
            if value=="application/json"
                for (k2, v2) in res.headers
                    if k2=="Content-Length" && v2=="0"
                        return HTTP.nobody
                    end
                end
                result = JSON.parse(read(IOBuffer(res.body), String);
                                                            dicttype=Dict{Symbol, Any})
                return raw || (res.status >= 400) ? result :
                                                method.transform(result, resource.transform)

            else
                result, status = res.body, res.status
                return status == 200 ? result : Dict{Symbol, Any}(:error =>
                                    Dict{Symbol, Any}(:message => result, :code => status))
            end
        end
    end
    nothing
end

const _default_session = Dict{APIRoot, GoogleSession}()

include("iam.jl")
include("storage.jl")
include("compute.jl")
include("container.jl")
include("pubsub.jl")
include("logging.jl")
include("datastore.jl")

end


text:
""
function constraint_power_balance_shunt_dispatch(pm::AbstractACPModel, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_gens, bus_storage, bus_shunts_var, bus_pd, bus_qd, bus_gs_const, bus_bs_const)
    vm   = var(pm, n, :vm, i)
    p    = get(var(pm, n),    :p, Dict()); PowerModels._check_var_keys(p, bus_arcs, "active power", "branch")
    q    = get(var(pm, n),    :q, Dict()); PowerModels._check_var_keys(q, bus_arcs, "reactive power", "branch")
    pg   = get(var(pm, n),   :pg, Dict()); PowerModels._check_var_keys(pg, bus_gens, "active power", "generator")
    qg   = get(var(pm, n),   :qg, Dict()); PowerModels._check_var_keys(qg, bus_gens, "reactive power", "generator")
    ps   = get(var(pm, n),   :ps, Dict()); PowerModels._check_var_keys(ps, bus_storage, "active power", "storage")
    qs   = get(var(pm, n),   :qs, Dict()); PowerModels._check_var_keys(qs, bus_storage, "reactive power", "storage")
    psw  = get(var(pm, n),  :psw, Dict()); PowerModels._check_var_keys(psw, bus_arcs_sw, "active power", "switch")
    qsw  = get(var(pm, n),  :qsw, Dict()); PowerModels._check_var_keys(qsw, bus_arcs_sw, "reactive power", "switch")
    p_dc = get(var(pm, n), :p_dc, Dict()); PowerModels._check_var_keys(p_dc, bus_arcs_dc, "active power", "dcline")
    q_dc = get(var(pm, n), :q_dc, Dict()); PowerModels._check_var_keys(q_dc, bus_arcs_dc, "reactive power", "dcline")

    bs = get(var(pm, n), :bs, Dict()); PowerModels._check_var_keys(bs, bus_shunts_var, "reactive power", "shunt")

    cstr_p = JuMP.@NLconstraint(pm.model, 0 == - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*vm^2)
    cstr_q = JuMP.@NLconstraint(pm.model, 0 == - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*vm^2 + sum(bs[s]*vm^2 for s in bus_shunts_var))

    if report_duals(pm)
        sol(pm, n, :bus, i)[:lam_kcl_r] = cstr_p
        sol(pm, n, :bus, i)[:lam_kcl_i] = cstr_q
    end
end

""
function constraint_power_balance_shunt_dispatch_soft(pm::AbstractACPModel, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_storage, bus_gens, bus_shunts_var, bus_pd, bus_qd, bus_gs_const, bus_bs_const)
    vm   = var(pm, n, :vm, i)
    p    = get(var(pm, n),    :p, Dict()); PowerModels._check_var_keys(p, bus_arcs, "active power", "branch")
    q    = get(var(pm, n),    :q, Dict()); PowerModels._check_var_keys(q, bus_arcs, "reactive power", "branch")
    pg   = get(var(pm, n),   :pg, Dict()); PowerModels._check_var_keys(pg, bus_gens, "active power", "generator")
    qg   = get(var(pm, n),   :qg, Dict()); PowerModels._check_var_keys(qg, bus_gens, "reactive power", "generator")
    ps   = get(var(pm, n),   :ps, Dict()); PowerModels._check_var_keys(ps, bus_storage, "active power", "storage")
    qs   = get(var(pm, n),   :qs, Dict()); PowerModels._check_var_keys(qs, bus_storage, "reactive power", "storage")
    psw  = get(var(pm, n),  :psw, Dict()); PowerModels._check_var_keys(psw, bus_arcs_sw, "active power", "switch")
    qsw  = get(var(pm, n),  :qsw, Dict()); PowerModels._check_var_keys(qsw, bus_arcs_sw, "reactive power", "switch")
    p_dc = get(var(pm, n), :p_dc, Dict()); PowerModels._check_var_keys(p_dc, bus_arcs_dc, "active power", "dcline")
    q_dc = get(var(pm, n), :q_dc, Dict()); PowerModels._check_var_keys(q_dc, bus_arcs_dc, "reactive power", "dcline")

    bs = get(var(pm, n), :bs, Dict()); PowerModels._check_var_keys(bs, bus_shunts_var, "reactive power", "shunt")

    p_delta_abs = var(pm, n, :p_delta_abs, i)
    q_delta_abs = var(pm, n, :q_delta_abs, i)

    p_delta = @NLexpression(pm.model, - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*vm^2)
    q_delta = @NLexpression(pm.model, - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*vm^2 + sum(bs[s]*vm^2 for s in bus_shunts_var))

    @NLconstraint(pm.model,  p_delta_abs >= p_delta)
    @NLconstraint(pm.model, -p_delta_abs <= p_delta)

    @NLconstraint(pm.model,  q_delta_abs >= q_delta)
    @NLconstraint(pm.model, -q_delta_abs <= q_delta)
end

""
function constraint_ohms_yt_from_goc(pm::AbstractACPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    p_fr  = var(pm, n,  :p, f_idx)
    q_fr  = var(pm, n,  :q, f_idx)
    vm_fr = var(pm, n, :vm, f_bus)
    vm_to = var(pm, n, :vm, t_bus)
    va_fr = var(pm, n, :va, f_bus)
    va_to = var(pm, n, :va, t_bus)

    JuMP.@NLconstraint(pm.model, p_fr ==  (g/tm^2+g_fr)*vm_fr^2 + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
    JuMP.@NLconstraint(pm.model, q_fr == -(b/tm^2+b_fr)*vm_fr^2 - (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
end

""
function expression_bus_withdrawal(pm::AbstractACPModel, n::Int, i::Int, bus_storage, bus_pd, bus_qd, bus_gs, bus_bs)
    vm = var(pm, n, :vm, i)
    ps = get(var(pm, n), :ps, Dict()); PowerModels._check_var_keys(ps, bus_storage, "active power", "storage")
    qs = get(var(pm, n), :qs, Dict()); PowerModels._check_var_keys(ps, bus_storage, "reactive power", "storage")

    ps_total = 0.0
    qs_total = 0.0
    if length(bus_storage) > 0
        ps_total = sum(ps[s] for s in bus_storage)
        qs_total = sum(qs[s] for s in bus_storage)
    end

    pd_total = 0.0
    qd_total = 0.0
    if length(bus_pd) > 0
        pd_total = sum(pd for pd in values(bus_pd))
    end
        if length(bus_qd) > 0
        qd_total = sum(qd for qd in values(bus_qd))
    end

    gs_total = 0.0
    bs_total = 0.0
    if length(bus_gs) > 0
        gs_total = sum(gs for gs in values(bus_gs))*vm^2
    end
    if length(bus_bs) > 0
        bs_total = sum(bs for bs in values(bus_bs))*vm^2
    end

    var(pm, n, :bus_wdp)[i] = ps_total + pd_total + gs_total
    var(pm, n, :bus_wdq)[i] = qs_total + qd_total + bs_total
end

""
function expression_branch_power_yt_from_goc(pm::AbstractACPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    vm_fr = var(pm, n, :vm, f_bus)
    vm_to = var(pm, n, :vm, t_bus)
    va_fr = var(pm, n, :va, f_bus)
    va_to = var(pm, n, :va, t_bus)

    var(pm, n, :p)[f_idx] = @NLexpression(pm.model,  (g/tm^2+g_fr)*vm_fr^2 + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
    var(pm, n, :q)[f_idx] = @NLexpression(pm.model, -(b/tm^2+b_fr)*vm_fr^2 - (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
end

""
function expression_branch_power_yt_from(pm::AbstractACPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    vm_fr = var(pm, n, :vm, f_bus)
    vm_to = var(pm, n, :vm, t_bus)
    va_fr = var(pm, n, :va, f_bus)
    va_to = var(pm, n, :va, t_bus)

    var(pm, n, :p)[f_idx] = @NLexpression(pm.model,  (g+g_fr)/tm^2*vm_fr^2 + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
    var(pm, n, :q)[f_idx] = @NLexpression(pm.model, -(b+b_fr)/tm^2*vm_fr^2 - (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
end

""
function expression_branch_power_yt_to(pm::AbstractACPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_to, b_to, tr, ti, tm)
    vm_fr = var(pm, n, :vm, f_bus)
    vm_to = var(pm, n, :vm, t_bus)
    va_fr = var(pm, n, :va, f_bus)
    va_to = var(pm, n, :va, t_bus)

    var(pm, n, :p)[t_idx] = @NLexpression(pm.model,  (g+g_to)*vm_to^2 + (-g*tr-b*ti)/tm^2*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/tm^2*(vm_to*vm_fr*sin(va_to-va_fr)) )
    var(pm, n, :q)[t_idx] = @NLexpression(pm.model, -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/tm^2*(vm_to*vm_fr*cos(va_to-va_fr)) + (-g*tr-b*ti)/tm^2*(vm_to*vm_fr*sin(va_to-va_fr)) )
end


text:
module ModelManager

export Model

using CoreExtended
using FileManager
using JLGEngine

using ..GraphicsManager
using ..StorageManager
using ..MeshManager

const Mesh = MeshManager.Typ

""" TODO """
type Model <: JLGEngine.IComponent
	id::Symbol
	name::String

	enabled::Bool
	linked::Bool
	loaded ::Bool

	source::FileSource

	meshes::Dict{Symbol, Mesh}
	group::StorageGroup
	mainStorage::StorageGroup

	Model(id=:_) = (g=StorageGroup(); new(id, "",true,false,false,FileSource(),Dict(),g,g))
end

""" TODO """
setMainBuffer(this::Model, storage::StorageGroup) = this.mainStorage = storage

""" TODO """
getMainBuffer(this::Model) = this.mainStorage

""" TODO """
link(this::Model) = StorageManager.link(this.group)

""" TODO """
function init(this::Model)
end

""" TODO """
function upload(this::Model)
	if this.loaded return end
	prepare(this::Model,this.meshes[:DEFAULT])
	StorageManager.upload(this.group)
	this.loaded=true
end

""" TODO """
function unload(this::Model)
	if !this.loaded return end
	# ....
	this.loaded=false
end

""" TODO """
function resets()
	GRAPHICSDRIVER().resetBuffers()
end

# preset for manager: includes code here
presetManager(Model)

#function create(k::Symbol, typ::Any)
#	if !haskey(list,k)
#		e=setData(typ)
#		list[k]=e
#		return e
#	end
#	list[k]
#end

""" TODO """
function getGroup(this::Model, k::Symbol)
	r = nothing
	if k == :VERTEX || k == :UV || k == :NORMAL
		r = StorageManager.getGroup(this.group.list[:VB].list[:VBO],k)
	elseif  k == :INDEX
		r = StorageManager.getGroup(this.group.list[:IB].list[:IBO],k)
	end
	r
	#StorageManager.getData(getMainBuffer(this),k)
end

""" TODO """
function getData(this::Model, k::Symbol)
	r = getGroup(this,k)
	if r != nothing r=r.data end #StorageManager.getValues
	r
end

""" TODO """
function getBuffer(this::Model, k::Symbol)
	#.....
end

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

include("ModelManager/MeshFabric.jl")

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

end #ModelManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
using Test
using IsURL

@testset "isurl" begin
    @test isurl("http://julialang.org") == true
    @test isurl("https://julialang.org") == true
    @test isurl("httpS://julialang.org") == true
    @test isurl("file://julialang.org") == true
    @test isurl("mailto:someone@example.com") == true
    @test isurl("data:text/plainbase64,SGVsbG8kIQ%3D") == true

    @test isurl("//julialang.org") == false
    @test isurl("/foo/bar") == false
    @test isurl("foo/bar") == false
    @test isurl("C:\\") == false
    @test isurl("c:\\Folder\\file") == false
    @test isurl("ht,tp://julialang.org") == false
end

@testset "isrelativeurl" begin
    @test isrelativeurl("../path/to/directory") == true
    @test isrelativeurl("./__file__") == true

    @test isrelativeurl("foo:bar") == false
    @test isrelativeurl("file://julialang.org") == false
end


text:
module Jacobi

export am,
    sn, cn, dn, nn,
    sd, cd, dd, nd,
    sc, cc, dc, nc,
    ss, cs, ds, ns

# Abramowitz & Stegun, section 16.4, p571
const _ambuf = Array{Float64}(undef, 10)
function _am(u::Float64, m::Float64, tol::Float64)
    if u == 0. return 0. end

    sqrt_tol = sqrt(tol)
    if m < sqrt_tol
        # A&S 16.13.4
        return u - 0.25*m*(u - 0.5*sin(2.0*u))
    end
    m1 = 1. - m
    if m1 < sqrt_tol
        # A&S 16.15.4
        t = tanh(u)
        return asin(t) + 0.25*m1*(t - u*(1. - t^2))*cosh(u)
    end

    a,b,c,n = 1., sqrt(m1), sqrt(m), 0
    while abs(c) > tol
        @assert n < 10
        a,b,c,n = 0.5*(a+b), sqrt(a*b), 0.5*(a-b), n+1
        _ambuf[n] = c/a
    end

    phi = ldexp(a*u, n)
    for i = n:-1:1
        phi = 0.5*(phi + asin(_ambuf[i]*sin(phi)))
    end
    phi
end
_am(u::Float64, m::Float64) = _am(u, m, eps(Float64))

"""
    am(u::Real, m::Real, [tol::Real=eps(Float64)])

Returns amplitude, φ, such that u = F(φ | m)

Landen sequence with convergence to `tol` used if `√(tol) ≤ m ≤ 1 - √(tol)`
"""
function am(u::Float64, m::Float64, tol::Float64)
    (m < 0. || m > 1.) && throw(DomainError(m, "argument m not in [0,1]"))
    _am(u, m, tol)
end
am(u::Float64, m::Float64) = am(u, m, eps(Float64))
am(u::Real, m::Real) = am(Float64(u), Float64(m))

for (f,a,b,c) in ((:sn, :(sin(phi)),                :(sqrtmu1*s), :(sqrt(mu)*sin(phi))),
                  (:cn, :(cos(phi)),                :(cos(phi)),  :(sqrt(1. - mu*sin(phi)^2))),
                  (:dn, :(sqrt(1. - m*sin(phi)^2)), :(1.),        :(cos(phi))))
    @eval begin
        function ($f)(u::Float64, m::Float64)
            # Abramowitz & Stegun, section 16.10, p573
            lt0 = m < 0.
            gt1 = m > 1.
            if !(lt0 || gt1)
                phi = _am(u,m)
                return $a
            elseif lt0
                mu1 = 1.0/(1. - m)
                mu = -m*mu1
                sqrtmu1 = sqrt(mu1)
                v = u/sqrtmu1
                phi = _am(v,mu)
                s = sin(phi)
                return ($b)/sqrt(1. - mu*s^2)
            elseif gt1
                mu = 1/m
                v = u*sqrt(m)
                phi = _am(v,mu)
                return $c
            end
        end
    end
end

xn = ((:s,:(sn(u,m))), (:c,:(cn(u,m))), (:d,:(dn(u,m))), (:n,:(1.)))
for (p,num) in xn, (q,den) in xn
    f = Symbol(p, q)
    @eval begin
        """
            $($f)(u::Real, m::Real)

        Compute the Jacobi elliptic function $($f)(u | m)
        """
        ($f)(u::Real, m::Real) = ($f)(Float64(u), Float64(m))
    end

    if (p == q)
        @eval ($f)(::Float64, ::Float64) = 1.0
    elseif (q != :n)
        @eval ($f)(u::Float64, m::Float64) = ($num)/($den)
    end
end

end # module


text:
using HMMs
using Brobdingnag
using Profile
using ProfileView
using GCTools

hmm1 = HMMs.randomhmm( HMMs.fullyconnected(2), calc=Brob )
y = HMMs.rand( hmm1, 100000 );

hmm2 = HMMs.randomhmm( HMMs.fullyconnected(2), calc=Brob, seed=2 )
HMMs.setobservations( hmm2, y );
@time HMMs.em( hmm2 )
hmm2 = HMMs.randomhmm( HMMs.fullyconnected(2), calc=Brob, seed=2 )
HMMs.setobservations( hmm2, y );
GCTools.reset()
@time HMMs.em( hmm2 )
GCTools.print()

Profile.clear()
hmm2 = HMMs.randomhmm( HMMs.fullyconnected(2), calc=Brob, seed=2 )
HMMs.setobservations( hmm2, y );
@profile HMMs.em( hmm2 )
ProfileView.view()

hmm2 = HMMs.randomhmm( HMMs.fullyconnected(2), calc=Brob, seed=2 )
HMMs.setobservations( hmm2, y );
@time HMMs.em( hmm2 )
# 2.02s
# 1.40s


text:
#%%
println("hello")
println("Julia")
#%%

#%%
num1=22
num2=22.35
println("n1: ",num1," ",typeof(num1))
println("n2: ",num2," ",typeof(num2))
d=num1/num2
println(d)
strA="myString"
println(typeof(strA))

num3=convert(Float64,num1)
println("n3: ",num3," ",typeof(num1))

#%%

#%%
function funcA(x)
    return x*5
end
println(funcA(2))
println(funcA([2,22]))
#%%