text:
## For parsing JSON into dlmreadable files
epitope_positions = [];
epitope_perauthor = Dict();
for f in readdir("./")
       if !isnothing(match(r".json", f))
              a = JSON.Parser.parsefile(f)
              epitope_perauthor[a["name"]] = []
              for (n,mp) in a["map"]
                     i0 = div(FluTools.gene_positions["h3n2","ha"][n][1], 3) + 1
                     for i in keys(mp)
                            push!(epitope_positions, i0+parse(Int, i))
                            push!(epitope_perauthor[a["name"]], i0+parse(Int,i))
                     end
                     println("$(a["name"]) / $(n) : $(length(mp)) positions")
              end
       end
end
       unique!(epitope_positions)
       sort!(epitope_positions)
       for (n,v) in epitope_perauthor
       sort!(v)
       end

## Loop for trajectories
param = get_aa_parameters()
authors = ["all", "koel", "luksza", "shih", "wolf"]
files = ["h3n2/data/epitopes/ha/aligned_h3n2_ha_aa_$(a)epitopes.fasta" for a in authors]
at_epitopes = Dict()
outplot_epitopes = Dict()
for (a,f) in zip(authors, files)
       println(a)
       at_epitopes[a] = get_trajectories(f, param)
       outplot_epitopes[a] = Pfix_vs_frequency(at_epitopes[a], param)
end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
module CJKFrequencies

using LightXML
using DataStructures
using LazyArtifacts

export CJKFrequency, charfreq,

SimplifiedLCMC, SimplifiedJunDa,

Lexicon, tagged_with

include("charfreq.jl")
include("frequency_datasets.jl")
include("lexicon.jl")

end

ext: jl


text:
function timings_check!(df::DataFrame; N = 6, upper=0.25)

    for row in 1:nrow(df)
        if df[row, :].max > (1.0 + upper) * df[row, :].median
            println("\n\nFixing row=$row: $(df[row, :])")

            #
            dft = timings(
                [df[row, :].num_threads], [df[row, :].num_cpp_chains],
                [df[row, :].num_chains], [df[row, :].num_samples], N)
            #
            df[row, :] = dft[1, :]
            println("Updating row=$row: $(df[row, :])")
        end
    end
end

#timings_check!(df)
#df

# Write df to results subdir

#CSV.write(joinpath(ProjDir, "results", "arm_no_tbb.csv"), df)
#CSV.write(joinpath(ProjDir, "results", "arm_tbb.csv"), df)
#CSV.write(joinpath(ProjDir, "results", "intel_no_tbb.csv"), df)
#CSV.write(joinpath(ProjDir, "results", "intel_tbb.csv"), df)

# Sort on (:num_threads,:num_chains0 instead on (:num_threads, :num_cpp_chains)

#sort(df, [:num_threads, :num_chains])

# Fix rows in df

function fix_row!(df::DataFrame, row)
    println(first(df, row))
    df[row, :] = timings(
        [df[row,:].num_threads], [df[row,:].num_cpp_chains],
        [df[row,:].num_chains], [df[row,:].num_samples], 6)[1, :]
    println("Updated row $row: $(df[row, :])")
end

#=
for i in 1:3
    fix_row!(df, i)
end
=#


text:
"""
`get_data`

Functions to find the path to data files according to character pattern considering that each session is a subfolder
Get data revised: this version operate in 2 possible way,
in the first way it will looks for file names containing a specified string
in the second way use a symbol to refer to a dictionary and find the specified string
"""
#Method 1 inputs a directory and a string
function get_data(dirnames,what::String)
    location = String[] #array to load with all the paths corrisponding to researched file type
    if eltype(dirnames) == Char #in case only one folder is loaded the for loop would research in single character
        tool = String[]
        push!(tool,dirnames)
        dirnames = tool
    end
    for dirname in dirnames
        files = readdir(dirname)
        for file in files
            if occursin(Regex(what), file) && !occursin(Regex(".txt"), file)
                complete_filename = joinpath(dirname,file)
                push!(location,complete_filename)
            end
        end
    end
    return location
end
#Method 2 a directory and a symbol
function get_data(dirnames, kind::Symbol)
    #the dictionary refers the symbol in the input to a specific string to look for
    ext_dict = Dict(:bhv => "a.csv", :cam => ".mat", :log => "AI.csv")
    if !(kind in keys(ext_dict))
        error("Choose among $(keys(ext_dict))")
    end
    #once the string is identified the function call itself again with the first method
    return get_data(dirnames, ext_dict[kind])
end

function get_data(Dir::String)
    files = readdir(Dir)
    bhv = []
    for file in files
        if occursin(Regex(".csv"), file)
            complete_filename = joinpath(Dir,file)
            push!(bhv,complete_filename)
        end
    end
    return bhv
end

"""
`create_DataIndex`
create a Dataframe to identify the raw files to processed it has 2 methods, find files in run_task_photo/raw_data
or takes all the files in a folder
"""
function create_DataIndex(bhv)
    string_search = match.(r"[a-zA-Z]{2}\d+_\d{6}[a-z]{1}",bhv)
    mask = string_search.!= nothing #remove files where there wasn't a match
    string_search = string_search[mask]
    bhv = bhv[mask]
    string_result = [res.match for res in string_search if res !== nothing]
    DataIndex = DataFrame(Bhv_Path = bhv)
    DataIndex[!,:Session] = String.(string_result.*".csv")
    DataIndex[!,:MouseID] = String.([split(t,"_")[1] for t in DataIndex[:,:Session]])
    DataIndex[!,:Day] = String.(["20"*match.(r"\d{6}",t).match for t in DataIndex[:,:Session]])
    DataIndex[!,:Period] = String.([match.(r"[a-z]{1}",t).match for t in DataIndex[:,:Session]])
    return DataIndex
end

"""
`find_behavior`
Function that deals with the type of preprocessing Single exp folder or raw data folder

"""

function find_behavior(Directory_path)
    Dir = replace(Directory_path,basename(Directory_path)=>"")
    saving_path = joinpath(Directory_path,"Bhv")
    if !ispath(saving_path)
        mkdir(saving_path)
    end
    bhv = get_data(Directory_path);
    DataIndex = create_DataIndex(bhv);
    DataIndex[!,:Preprocessed_Path] = joinpath.(saving_path,DataIndex[:,:Session])
    results_path = joinpath(Directory_path,"Results")
    if !ispath(results_path)
        mkdir(results_path)
    end
    DataIndex[!,:Saving_Path] .= results_path
    return DataIndex
end

function find_behavior(Directory_path::String, Exp_type::String,Exp_name::String, Mice_suffix ::String; run_task = "run_task_photo")
    rawdata_path = joinpath(Directory_path,run_task,"raw_data")
    exp_dir = joinpath(Directory_path,"Datasets",Exp_type,Exp_name)
    saving_path = joinpath(Directory_path,"Datasets",Exp_type,Exp_name,"Bhv")
    if !ispath(saving_path)
        if !ispath(exp_dir)
            mkdir(exp_dir)
        end
        mkdir(saving_path)
    end
    bhv = get_data(rawdata_path, Mice_suffix);
    DataIndex = create_DataIndex(bhv);
    DataIndex[!,:Saving_Path] .= exp_dir
    DataIndex[!,:Preprocessed_Path] = joinpath.(saving_path,DataIndex[:,:Session])
    return DataIndex
end

"""
`paths_dataframe`
Create a Dataframe to store paths of files to preprocess
"""
function paths_dataframe(bhv)
    behavior = DataFrame()
    mask = occursin.(bhv,"txt")
    bhv = bhv[mask]
    ##### extract date and mouse ID per session using get_mousedate (it works with a full path)
    MouseID = Vector{String}(undef, size(bhv,1))
    Day = Vector{String}(undef, size(bhv,1))
    Session = Vector{String}(undef, size(bhv,1))
    for i = 1:size(bhv,1)
        MouseID[i], Day[i], Session[i] = Flipping.get_BHVmousedate(bhv[i])
    end
    behavior = DataFrame(MouseID = MouseID)
    behavior[:Day] = Day #file properties are not reliable for the date of the session
    behavior[:Session] = Session.*".csv"
    behavior[:Bhv_Path] = bhv

    return behavior
end

"""
`createfilelist`
use get_data function to obtain all filenames of behaviour
"""
function createfilelist(Directory_path::String, Mice_suffix::String)
    bhv = get_data(Directory_path,:bhv)
    #use get_sessionname to select relevant session (for instance use exp naming code)
    bhv_session = map(t->get_sessionname(t,Mice_suffix),bhv)#to be changed for each dataframe
    # get_sessionname return a start result for sessions that don't match the criteria this can be used to prune irrelevant paths
    bhv = bhv[bhv_session.!="start"]
    bhv_session = bhv_session[bhv_session.!="start"]
    return bhv
end

"""
`get_sessionname(filepath, what::String)`
Use it to find the name of a session from a path, can accept a string or a symbol connected to a dict to find
file matching the requirements
"""
#this function extract the name of a session from a filepath according to a given character pattern
function get_sessionname(filepath, what::String)
    pathinfo = split(filepath,"/")
    sessionname = "start"
    for piece in pathinfo
        if occursin(Regex(what), string(piece))
            sessionname = piece
        end
    end
    sessionname = String(sessionname)
    return sessionname
end

# the second method allows to save experiment related character pattern in a dictionary
function get_sessionname(filepath, kind::Symbol)
    #the dictionary refers the symbol in the input to a specific string to look for
    ext_dict = Dict(:GcAMP => "170", :BilNac => "NB")
    if !(kind in keys(ext_dict))
        error("Choose among $(keys(ext_dict))")
    end
    #once the string is identified the function call itself again with the first method
    return get_sessionname(filepath, ext_dict[kind])
end

"""
`get_session`
Generalise version
"""
function get_session(filepath,what::String)
    a = get_sessionname(filepath,what)
    b = match(r"[a-zA-Z]{2}\d+",a)
    c = b.match
    if occursin(r"\d{8}",a)
        d = match(r"\d{8}",a)
        e = d.match
        e = e[3:8]
    else
        d = match(r"\d{6}",a)
        e = c.match
    end
    sessione = c*"_"*e
    return sessione
end

"""
`check_fiberlocation`

look for a dataset where fiberlocation across day is stored
"""

function check_fiberlocation(data,exp_dir)
    filetofind=joinpath(exp_dir,"FiberLocation.csv");
    if isfile(filetofind)
        fiberlocation = FileIO.load(filetofind) |>DataFrame;
        merged_data = join(data, fiberlocation, on = :Session, kind = :left,makeunique=true);
        println("found fibres location file, HAVE YOU UPDATED IT?")
        return merged_data
    else
        println("no fibres location file")
        return data
    end
end

function check_fiberlocation(data,Directory_path,Exp_name;run_task = "run_task_photo")
    exp_dir = joinpath(Directory_path,run_task,Exp_name)
    check_fiberlocation(data,exp_dir)
end


text:
using PyPlot
using JLD

include("model_tyc3559_v2.jl")
include("transit_orb.jl")
include("occultquad.jl")
include("kepler.jl")
include("compute_likelihood_save.jl")
include("lorentz_likelihood_hermitian_band_save.jl")
include("banbks.jl")

@load "tyc3559_mcmc_norv_fixedkernel_v02.jld"

# Initialize walkers from last step of prior run:
for j=1:nwalkers
  log_like_mcmc[j,1]=log_like_mcmc[j,nsteps]
  par_mcmc[j,1,:]=par_mcmc[j,nsteps,:]
end

# Initialize scale length & acceptance counter:
ascale = 2.0
accept = 0
tic()
# Next, loop over steps in markov chain:
for i=2:nsteps
  for j=1:nwalkers
    ipartner = j
# Choose another walker to 'breed' a trial step with:
    while ipartner == j
      ipartner = ceil(Int,rand()*nwalkers)
    end
# Now, choose a new set of parameters for the trial step:
    z=(rand()*(sqrt(ascale)-1.0/sqrt(ascale))+1.0/sqrt(ascale))^2
    par_trial = vec(par_mcmc[j,i-1,:])
    for itrial=1:nvary
      par_trial[ivary[itrial]]=z*par_mcmc[j,i-1,ivary[itrial]]+(1.0-z)*par_mcmc[ipartner,i-1,ivary[itrial]]
    end
# Compute model & chi-square:
    aonr2 = (1./par_trial[5]*p0/pi)^2*(1-par_trial[3]^2)
    if (par_trial[6] < 0) || (par_trial[6] > 1) || (par_trial[7] < 0) || (par_trial[7] > 1) || (par_trial[4] < 0) || (par_trial[5] < 0) || (abs(par_trial[3]) > 1) || (par_trial[3]^2 > aonr2) then
      log_like_trial = -1e100
    else
      log_like_trial = compute_likelihood_save(time_phot,par_trial,flux,aex,al_small,indx,logdeta,p)
    end
# Next, determine whether to accept this trial step:
    alp = z^(nvary-1)*exp(log_like_trial - log_like_mcmc[j,i-1])
    if alp >= rand()
# If step is accepted, add it to the chains!
      par_mcmc[j,i,:] = par_trial
      log_like_mcmc[j,i,:] = log_like_trial
      accept = accept + 1
    else
# If step is rejected, then copy last step:
      par_mcmc[j,i,:] = par_mcmc[j,i-1,:]
      log_like_mcmc[j,i,:] = log_like_mcmc[j,i-1]
    end
  end
  if mod(i,100) == 0
    eltime = toq()
    println("Number of steps: ",i," acceptance rate: ",accept/(100*nwalkers)," elapsed time [hr]: ",eltime/3600.)
    accept = 0
    tic()
  end
end

for i=1:nvary
  for j=1:nwalkers
    plot(vec(par_mcmc[j,1:nsteps,ivary[i]]))
  end
  xlabel("MCMC step")
  ylabel(pname[ivary[i]])
  println("Hit return to continue")
  read(STDIN,Char)
  clf()
end

@save "tyc3559_mcmc_norv_fixedkernel_v03.jld"


text:
const _DESERIALIZE_MESSAGE = "Deserialized initial_conditions_data"
const _MAKE_IC_MESSAGE = "Make Initial Conditions Model"
const _SKIP_IC_MESSAGE = "Skip build of initial conditions"

function test_ic_serialization_outputs(model::PSI.OperationModel; ic_file_exists, message)
    ic_file = PSI.get_initial_conditions_file(model)
    log_file = PSI.get_log_file(model)

    @test isfile(ic_file) == ic_file_exists
    if ic_file_exists
        @test Serialization.deserialize(ic_file) isa PSI.InitialConditionsData
    end

    make = false
    deserialize = false
    skip = false
    if message == "make"
        make = true
    elseif message == "deserialize"
        deserialize = true
    elseif message == "skip"
        skip = true
    else
        error("invalid: $message")
    end

    text = read(log_file, String)
    @test make == occursin(_MAKE_IC_MESSAGE, text)
    @test deserialize == occursin(_DESERIALIZE_MESSAGE, text)
    @test skip == occursin(_SKIP_IC_MESSAGE, text)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ReplicationControllerSpec is the specification of a replication controller.

    IoK8sApiCoreV1ReplicationControllerSpec(;
        minReadySeconds=nothing,
        replicas=nothing,
        selector=nothing,
        template=nothing,
    )

    - minReadySeconds::Int32 : Minimum number of seconds for which a newly created pod should be ready without any of its container crashing, for it to be considered available. Defaults to 0 (pod will be considered available as soon as it is ready)
    - replicas::Int32 : Replicas is the number of desired replicas. This is a pointer to distinguish between explicit zero and unspecified. Defaults to 1. More info: https://kubernetes.io/docs/concepts/workloads/controllers/replicationcontroller#what-is-a-replicationcontroller
    - selector::Dict{String, String} : Selector is a label query over pods that should match the Replicas count. If Selector is empty, it is defaulted to the labels present on the Pod template. Label keys and values that must match in order to be controlled by this replication controller, if empty defaulted to labels on Pod template. More info: https://kubernetes.io/docs/concepts/overview/working-with-objects/labels/#label-selectors
    - template::IoK8sApiCoreV1PodTemplateSpec : Template is the object that describes the pod that will be created if insufficient replicas are detected. This takes precedence over a TemplateRef. More info: https://kubernetes.io/docs/concepts/workloads/controllers/replicationcontroller#pod-template
"""
mutable struct IoK8sApiCoreV1ReplicationControllerSpec <: SwaggerModel
    minReadySeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: minReadySeconds
    replicas::Any # spec type: Union{ Nothing, Int32 } # spec name: replicas
    selector::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: selector
    template::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodTemplateSpec } # spec name: template

    function IoK8sApiCoreV1ReplicationControllerSpec(;minReadySeconds=nothing, replicas=nothing, selector=nothing, template=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ReplicationControllerSpec, Symbol("minReadySeconds"), minReadySeconds)
        setfield!(o, Symbol("minReadySeconds"), minReadySeconds)
        validate_property(IoK8sApiCoreV1ReplicationControllerSpec, Symbol("replicas"), replicas)
        setfield!(o, Symbol("replicas"), replicas)
        validate_property(IoK8sApiCoreV1ReplicationControllerSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiCoreV1ReplicationControllerSpec, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        o
    end
end # type IoK8sApiCoreV1ReplicationControllerSpec

const _property_map_IoK8sApiCoreV1ReplicationControllerSpec = Dict{Symbol,Symbol}(Symbol("minReadySeconds")=>Symbol("minReadySeconds"), Symbol("replicas")=>Symbol("replicas"), Symbol("selector")=>Symbol("selector"), Symbol("template")=>Symbol("template"))
const _property_types_IoK8sApiCoreV1ReplicationControllerSpec = Dict{Symbol,String}(Symbol("minReadySeconds")=>"Int32", Symbol("replicas")=>"Int32", Symbol("selector")=>"Dict{String, String}", Symbol("template")=>"IoK8sApiCoreV1PodTemplateSpec")
Base.propertynames(::Type{ IoK8sApiCoreV1ReplicationControllerSpec }) = collect(keys(_property_map_IoK8sApiCoreV1ReplicationControllerSpec))
Swagger.property_type(::Type{ IoK8sApiCoreV1ReplicationControllerSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ReplicationControllerSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ReplicationControllerSpec }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ReplicationControllerSpec[property_name]

function check_required(o::IoK8sApiCoreV1ReplicationControllerSpec)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ReplicationControllerSpec }, name::Symbol, val)
end


text:
export
    issue008,
    issue010,
    issue100,
    issue143,
    issue146a,
    issue146b,
    issue176

function issue008()
    large = [i^4 + 1e12 for i in 1:10^3]
    small = [1:10^3;]
    plot(small,large)
end

function issue010()
    p = plot(abs.(sin.(0:.1:10)))
    setattr(p, "ylog", true)
    setattr(p, "yrange", (5e-6,0.2))
    p
end

function issue100()
    x = range(0, stop=6pi, length=100)
    y = sin.(x)
    stem(x, y, "b;")
end

function issue143()
    p = plot(1:10, yrange=(10,0))
    setattr(p.y2, draw_ticklabels=true)
    p
end

issue146a() = plot(sin, 0, 2pi)
issue146b() = plot(sin, cos, xrange=[0,2pi], yrange=[-2,2])

function issue176()
    t = Table(2,2)
    z = reshape(1:100, 10, 10)
    t[1,1] = imagesc((1,10), (1,10), z)
    t[1,2] = imagesc((1,10), (10,1), z)
    t[2,1] = imagesc((10,1), (1,10), z)
    t[2,2] = imagesc((10,1), (10,1), z)
    t
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiNetworkingV1NetworkPolicySpec <: SwaggerModel
    egress::Any # spec type: Union{ Nothing, Vector{IoK8sApiNetworkingV1NetworkPolicyEgressRule} } # spec name: egress
    ingress::Any # spec type: Union{ Nothing, Vector{IoK8sApiNetworkingV1NetworkPolicyIngressRule} } # spec name: ingress
    podSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: podSelector
    policyTypes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: policyTypes

    function IoK8sApiNetworkingV1NetworkPolicySpec(;egress=nothing, ingress=nothing, podSelector=nothing, policyTypes=nothing)
        o = new()
        validate_property(IoK8sApiNetworkingV1NetworkPolicySpec, Symbol("egress"), egress)
        setfield!(o, Symbol("egress"), egress)
        validate_property(IoK8sApiNetworkingV1NetworkPolicySpec, Symbol("ingress"), ingress)
        setfield!(o, Symbol("ingress"), ingress)
        validate_property(IoK8sApiNetworkingV1NetworkPolicySpec, Symbol("podSelector"), podSelector)
        setfield!(o, Symbol("podSelector"), podSelector)
        validate_property(IoK8sApiNetworkingV1NetworkPolicySpec, Symbol("policyTypes"), policyTypes)
        setfield!(o, Symbol("policyTypes"), policyTypes)
        o
    end
end # type IoK8sApiNetworkingV1NetworkPolicySpec

const _property_map_IoK8sApiNetworkingV1NetworkPolicySpec = Dict{Symbol,Symbol}(Symbol("egress")=>Symbol("egress"), Symbol("ingress")=>Symbol("ingress"), Symbol("podSelector")=>Symbol("podSelector"), Symbol("policyTypes")=>Symbol("policyTypes"))
const _property_types_IoK8sApiNetworkingV1NetworkPolicySpec = Dict{Symbol,String}(Symbol("egress")=>"Vector{IoK8sApiNetworkingV1NetworkPolicyEgressRule}", Symbol("ingress")=>"Vector{IoK8sApiNetworkingV1NetworkPolicyIngressRule}", Symbol("podSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("policyTypes")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiNetworkingV1NetworkPolicySpec }) = collect(keys(_property_map_IoK8sApiNetworkingV1NetworkPolicySpec))
Swagger.property_type(::Type{ IoK8sApiNetworkingV1NetworkPolicySpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiNetworkingV1NetworkPolicySpec[name]))}
Swagger.field_name(::Type{ IoK8sApiNetworkingV1NetworkPolicySpec }, property_name::Symbol) =  _property_map_IoK8sApiNetworkingV1NetworkPolicySpec[property_name]

function check_required(o::IoK8sApiNetworkingV1NetworkPolicySpec)
    (getproperty(o, Symbol("podSelector")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiNetworkingV1NetworkPolicySpec }, name::Symbol, val)
end


text:
module PEPOKit

using TensorKit
using LinearAlgebra
using KrylovKit
using MacroTools: @forward

import TensorKit: Group, AbelianGroup

export AbstractPEPO,
    PEPOTensor,
    FinitePEPO,
    InfinitePEPO,
    AbstractPhysicalTensor,
    Liouvillian,
    SymmetricLiouvillian,
    StandardLiouvillian
export VectorPEPOTensor
export bloch_pepo
export cj, conjugate

export σ0, σ1, σ2, σ3, σp, σm
export SymmetryGroup

include("physicaltensor.jl")
include("pepo.jl")
include("vectorpepo.jl")
include("linalg.jl")
include("lindblad.jl")
include("symmetrise.jl")
include("tensoroperations.jl")
include("fet/vumps.jl")
include("fet/bondenv.jl")
include("fet/isometries.jl")

end # module

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


text:
# This file is based on a part of Julia. License is MIT: https://julialang.org/license

using QuasiArrays, Test, LinearAlgebra
import QuasiArrays: MemoryLayout

@testset "QuasiAdjoint/Transpose" begin
    @testset "Adjoint and Transpose inner constructor basics" begin
        intvec, intmat = QuasiVector([1, 2],0:0.5:0.5), QuasiArray([1 2; 3 4],(0:0.5:0.5,0:0.5:0.5))
        # Adjoint/Transpose inner constructor basic functionality, concrete scalar eltype
        @test (QuasiAdjoint{Int,typeof(intvec)}(intvec)::QuasiAdjoint{Int,typeof(intvec)}).parent === intvec
        @test (QuasiAdjoint{Int,typeof(intmat)}(intmat)::QuasiAdjoint{Int,typeof(intmat)}).parent === intmat
        @test (QuasiTranspose{Int,typeof(intvec)}(intvec)::QuasiTranspose{Int,typeof(intvec)}).parent === intvec
        @test (QuasiTranspose{Int,typeof(intmat)}(intmat)::QuasiTranspose{Int,typeof(intmat)}).parent === intmat
        # Adjoint/Transpose inner constructor basic functionality, abstract scalar eltype
        anyvec, anymat = QuasiVector(Any[1, 2],0:0.5:0.5), QuasiArray(Any[1 2; 3 4],(0:0.5:0.5,0:0.5:0.5))
        @test (QuasiAdjoint{Any,typeof(anyvec)}(anyvec)::QuasiAdjoint{Any,typeof(anyvec)}).parent === anyvec
        @test (QuasiAdjoint{Any,typeof(anymat)}(anymat)::QuasiAdjoint{Any,typeof(anymat)}).parent === anymat
        @test (QuasiTranspose{Any,typeof(anyvec)}(anyvec)::QuasiTranspose{Any,typeof(anyvec)}).parent === anyvec
        @test (QuasiTranspose{Any,typeof(anymat)}(anymat)::QuasiTranspose{Any,typeof(anymat)}).parent === anymat
    end

    @testset "Adjoint and Transpose outer constructor basics" begin
        intvec, intmat = QuasiVector([1, 2],0:0.5:0.5), QuasiArray([1 2; 3 4],(0:0.5:0.5,0:0.5:0.5))
        # the wrapped array's eltype strictly determines the Adjoint/Transpose eltype
        # so Adjoint{T}/Transpose{T} constructors are somewhat unnecessary and error-prone
        # so ascertain that such calls throw whether or not T and the input eltype are compatible
        @test_throws MethodError QuasiAdjoint{Int}(intvec)
        @test_throws MethodError QuasiAdjoint{Int}(intmat)
        @test_throws MethodError QuasiAdjoint{Float64}(intvec)
        @test_throws MethodError QuasiAdjoint{Float64}(intmat)
        @test_throws MethodError QuasiTranspose{Int}(intvec)
        @test_throws MethodError QuasiTranspose{Int}(intmat)
        @test_throws MethodError QuasiTranspose{Float64}(intvec)
        @test_throws MethodError QuasiTranspose{Float64}(intmat)
        # QuasiAdjoint/QuasiTranspose outer constructor basic functionality, concrete scalar eltype
        @test (QuasiAdjoint(intvec)::QuasiAdjoint{Int,typeof(intvec)}).parent === intvec
        @test (QuasiAdjoint(intmat)::QuasiAdjoint{Int,typeof(intmat)}).parent === intmat
        @test (QuasiTranspose(intvec)::QuasiTranspose{Int,typeof(intvec)}).parent === intvec
        @test (QuasiTranspose(intmat)::QuasiTranspose{Int,typeof(intmat)}).parent === intmat
        # the tests for the inner constructors exercise abstract scalar and concrete array eltype, forgoing here
    end

    @testset "Adjoint and Transpose add additional layers to already-wrapped objects" begin
        intvec, intmat = QuasiVector([1, 2],0:0.5:0.5), QuasiArray([1 2; 3 4],(0:0.5:0.5,0:0.5:0.5))
        @test (A = QuasiAdjoint(QuasiAdjoint(intvec))::QuasiAdjoint{Int,QuasiAdjoint{Int,typeof(intvec)}}; A.parent.parent === intvec)
        @test (A = QuasiAdjoint(QuasiAdjoint(intmat))::QuasiAdjoint{Int,QuasiAdjoint{Int,typeof(intmat)}}; A.parent.parent === intmat)
        @test (A = Transpose(Transpose(intvec))::Transpose{Int,Transpose{Int,typeof(intvec)}}; A.parent.parent === intvec)
        @test (A = Transpose(Transpose(intmat))::Transpose{Int,Transpose{Int,typeof(intmat)}}; A.parent.parent === intmat)
    end

    @testset "Adjoint and Transpose basic AbstractArray functionality" begin
        # vectors and matrices with real scalar eltype, and their adjoints/transposes
        intvec, intmat = QuasiVector([1, 2],0:0.5:0.5), QuasiMatrix([1 2 3; 4 5 6],(0:0.5:0.5, Base.OneTo(3)))
        tintvec, tintmat = QuasiMatrix([1 2],(Base.OneTo(1),0:0.5:0.5)), QuasiMatrix([1 4; 2 5; 3 6],(Base.OneTo(3),0:0.5:0.5))
        @testset "length methods" begin
            @test length(QuasiAdjoint(intvec)) == length(intvec)
            @test length(QuasiAdjoint(intmat)) == length(intmat)
            @test length(QuasiTranspose(intvec)) == length(intvec)
            @test length(Transpose(intmat)) == length(intmat)
        end
        @testset "size methods" begin
            @test size(QuasiAdjoint(intvec)) == (1, length(intvec))
            @test size(QuasiAdjoint(intmat)) == reverse(size(intmat))
            @test size(QuasiTranspose(intvec)) == (1, length(intvec))
            @test size(QuasiTranspose(intmat)) == reverse(size(intmat))
        end
        @testset "indices methods" begin
            @test axes(QuasiAdjoint(intvec)) == (Base.OneTo(1), Inclusion(0:0.5:0.5))
            @test axes(QuasiAdjoint(intmat)) == reverse(axes(intmat))
            @test axes(QuasiTranspose(intvec)) == (Base.OneTo(1), Inclusion(0:0.5:0.5))
            @test axes(QuasiTranspose(intmat)) == reverse(axes(intmat))
        end
        @testset "IndexStyle methods" begin
            @test IndexStyle(QuasiAdjoint(intvec)) == IndexCartesian()
            @test IndexStyle(QuasiAdjoint(intmat)) == IndexCartesian()
            @test IndexStyle(QuasiTranspose(intvec)) == IndexCartesian()
            @test IndexStyle(QuasiTranspose(intmat)) == IndexCartesian()
        end
        # vectors and matrices with complex scalar eltype, and their adjoints/transposes
        complexintvec, complexintmat = QuasiVector([1im, 2im],0:0.5:0.5), QuasiMatrix([1im 2im 3im; 4im 5im 6im],(0:0.5:0.5, Base.OneTo(3)))
        tcomplexintvec, tcomplexintmat = QuasiMatrix([1im 2im],(Base.OneTo(1),0:0.5:0.5)), QuasiMatrix([1im 4im; 2im 5im; 3im 6im],(Base.OneTo(3),0:0.5:0.5))
        acomplexintvec, acomplexintmat = conj.(tcomplexintvec), conj.(tcomplexintmat)
        # vectors and matrices with real-vector and real-matrix eltype, and their adjoints/transposes
        intvecvec = QuasiVector([[1, 2], [3, 4]],0:0.5:0.5)
        tintvecvec = QuasiMatrix([[[1 2]] [[3 4]]], (Base.OneTo(1),0:0.5:0.5))
        intmatmat = QuasiMatrix([[[1 2]] [[3  4]] [[ 5  6]];
                    [[7 8]] [[9 10]] [[11 12]]], (0:0.5:0.5,Base.OneTo(3)))
        tintmatmat = QuasiMatrix([[hcat([1, 2])] [hcat([7, 8])];
                    [hcat([3, 4])] [hcat([9, 10])];
                    [hcat([5, 6])] [hcat([11, 12])]], (Base.OneTo(3),0:0.5:0.5))
        # vectors and matrices with complex-vector and complex-matrix eltype, and their adjoints/transposes
        complexintvecvec, complexintmatmat = im .* (intvecvec, intmatmat)
        tcomplexintvecvec, tcomplexintmatmat = im .* (tintvecvec, tintmatmat)
        acomplexintvecvec, acomplexintmatmat = conj.(tcomplexintvecvec), conj.(tcomplexintmatmat)
        @testset "getindex methods, elementary" begin
            # implicitly test elementary definitions, for arrays with concrete real scalar eltype
            @test axes(QuasiAdjoint(intvec)) == axes(tintvec)
            @test QuasiAdjoint(intvec) == tintvec
            @test QuasiAdjoint(intmat) == tintmat
            @test QuasiTranspose(intvec) == tintvec
            @test QuasiTranspose(intmat) == tintmat
            # implicitly test elementary definitions, for arrays with concrete complex scalar eltype
            @test QuasiAdjoint(complexintvec) == acomplexintvec
            @test QuasiAdjoint(complexintmat) == acomplexintmat
            @test QuasiTranspose(complexintvec) == tcomplexintvec
            @test QuasiTranspose(complexintmat) == tcomplexintmat
            # implicitly test elementary definitions, for arrays with concrete real-array eltype
            @test QuasiAdjoint(intvecvec) == tintvecvec
            @test QuasiAdjoint(intmatmat) == tintmatmat
            @test QuasiTranspose(intvecvec) == tintvecvec
            @test QuasiTranspose(intmatmat) == tintmatmat
            # implicitly test elementary definitions, for arrays with concrete complex-array type
            @test QuasiAdjoint(complexintvecvec) == acomplexintvecvec
            @test QuasiAdjoint(complexintmatmat) == acomplexintmatmat
            @test QuasiTranspose(complexintvecvec) == tcomplexintvecvec
            @test QuasiTranspose(complexintmatmat) == tcomplexintmatmat
        end
        @testset "getindex(::AdjOrTransVec, ::Colon, ::AbstractArray{Int}) methods that preserve wrapper type" begin
            # for arrays with concrete scalar eltype
            @test QuasiAdjoint(intvec)[:, [0.0, 0.5]] == Adjoint(parent(intvec))
            @test QuasiTranspose(intvec)[:, [0.0, 0.5]] == Transpose(parent(intvec))
            @test QuasiAdjoint(complexintvec)[:, [0.0, 0.5]] == Adjoint(parent(complexintvec))
            @test QuasiTranspose(complexintvec)[:, [0.0, 0.5]] == Transpose(parent(complexintvec))
            # for arrays with concrete array eltype
            @test QuasiAdjoint(intvecvec)[:, [0.0, 0.5]] == Adjoint(parent(intvecvec))
            @test QuasiTranspose(intvecvec)[:, [0.0, 0.5]] == Transpose(parent(intvecvec))
            @test QuasiAdjoint(complexintvecvec)[:, [0.0, 0.5]] == Adjoint(parent(complexintvecvec))
            @test QuasiTranspose(complexintvecvec)[:, [0.0, 0.5]] == Transpose(parent(complexintvecvec))
        end
        @testset "getindex(::AdjOrTransVec, ::Colon, ::Colon) methods that preserve wrapper type" begin
            # for arrays with concrete scalar eltype
            @test QuasiAdjoint(intvec)[:, :] == QuasiAdjoint(intvec)
            @test QuasiTranspose(intvec)[:, :] == QuasiTranspose(intvec)
            @test QuasiAdjoint(complexintvec)[:, :] == QuasiAdjoint(complexintvec)
            @test QuasiTranspose(complexintvec)[:, :] == QuasiTranspose(complexintvec)
            # for arrays with concrete array elype
            @test QuasiAdjoint(intvecvec)[:, :] == QuasiAdjoint(intvecvec)
            @test QuasiTranspose(intvecvec)[:, :] == QuasiTranspose(intvecvec)
            @test QuasiAdjoint(complexintvecvec)[:, :] == QuasiAdjoint(complexintvecvec)
            @test QuasiTranspose(complexintvecvec)[:, :] == QuasiTranspose(complexintvecvec)
        end
        @testset "getindex(::AdjOrTransVec, ::Colon, ::Int) should preserve wrapper type on result entries" begin
            # for arrays with concrete scalar eltype
            @test QuasiAdjoint(intvec)[:, 0.5] == intvec[0.5:0.5]
            @test QuasiTranspose(intvec)[:, 0.5] == intvec[0.5:0.5]
            @test QuasiAdjoint(complexintvec)[:, 0.5] == conj.(complexintvec[0.5:0.5])
            @test QuasiTranspose(complexintvec)[:, 0.5] == complexintvec[0.5:0.5]
            # for arrays with concrete array eltype
            @test QuasiAdjoint(intvecvec)[:, 0.5] == Adjoint.(intvecvec[0.5:0.5])
            @test QuasiTranspose(intvecvec)[:, 0.5] == Transpose.(intvecvec[0.5:0.5])
            @test QuasiAdjoint(complexintvecvec)[:, 0.5] == Adjoint.(complexintvecvec[0.5:0.5])
            @test QuasiTranspose(complexintvecvec)[:, 0.5] == Transpose.(complexintvecvec[0.5:0.5])
        end
        @testset "setindex! methods" begin
            # for vectors with real scalar eltype
            @test (wv = QuasiAdjoint(copy(intvec));
                    wv === setindex!(wv, 3, 1, 0.5) &&
                    wv == setindex!(copy(tintvec), 3, 1, 0.5)    )
            @test (wv = QuasiTranspose(copy(intvec));
                    wv === setindex!(wv, 4, 1, 0.5) &&
                    wv == setindex!(copy(tintvec), 4, 1, 0.5)    )
            # for matrices with real scalar eltype
            @test (wA = QuasiAdjoint(copy(intmat));
                    wA === setindex!(wA, 7, 3, 0.0) &&
                    wA == setindex!(copy(tintmat), 7, 3, 0.0)    )
            @test (wA = QuasiTranspose(copy(intmat));
                    wA === setindex!(wA, 7, 3, 0.0) &&
                    wA == setindex!(copy(tintmat), 7, 3, 0.0)    )
            # for vectors with complex scalar eltype
            @test (wz = QuasiAdjoint(copy(complexintvec));
                    wz === setindex!(wz, 3im, 1, 0.5) &&
                    wz == setindex!(copy(acomplexintvec), 3im, 1, 0.5)   )
            @test (wz = QuasiTranspose(copy(complexintvec));
                    wz === setindex!(wz, 4im, 1, 0.5) &&
                    wz == setindex!(copy(tcomplexintvec), 4im, 1, 0.5)   )
            # for  matrices with complex scalar eltype
            @test (wZ = QuasiAdjoint(copy(complexintmat));
                    wZ === setindex!(wZ, 7im, 3, 0.0) &&
                    wZ == setindex!(copy(acomplexintmat), 7im, 3, 0.0)   )
            @test (wZ = QuasiTranspose(copy(complexintmat));
                    wZ === setindex!(wZ, 7im, 3, 0.0) &&
                    wZ == setindex!(copy(tcomplexintmat), 7im, 3, 0.0)   )
            # for vectors with concrete real-vector eltype
            @test (wv = QuasiAdjoint(copy(intvecvec));
                    wv === setindex!(wv, Adjoint([5, 6]), 1, 0.5) &&
                    wv == setindex!(copy(tintvecvec), [5 6], 1, 0.5))
            @test (wv = QuasiTranspose(copy(intvecvec));
                    wv === setindex!(wv, Transpose([5, 6]), 1, 0.5) &&
                    wv == setindex!(copy(tintvecvec), [5 6], 1, 0.5))
            # for matrices with concrete real-matrix eltype
            @test (wA = QuasiAdjoint(copy(intmatmat));
                    wA === setindex!(wA, Adjoint([13 14]), 3, 0.5) &&
                    wA == setindex!(copy(tintmatmat), hcat([13, 14]), 3, 0.5))
            @test (wA = QuasiTranspose(copy(intmatmat));
                    wA === setindex!(wA, Transpose([13 14]), 3, 0.5) &&
                    wA == setindex!(copy(tintmatmat), hcat([13, 14]), 3, 0.5))
            # for vectors with concrete complex-vector eltype
            @test (wz = QuasiAdjoint(copy(complexintvecvec));
                    wz === setindex!(wz, Adjoint([5im, 6im]), 1, 0.5) &&
                    wz == setindex!(copy(acomplexintvecvec), [-5im -6im], 1, 0.5))
            @test (wz = QuasiTranspose(copy(complexintvecvec));
                    wz === setindex!(wz, QuasiTranspose([5im, 6im]), 1, 0.5) &&
                    wz == setindex!(copy(tcomplexintvecvec), [5im 6im], 1, 0.5))
            # for matrices with concrete complex-matrix eltype
            @test (wZ = QuasiAdjoint(copy(complexintmatmat));
                    wZ === setindex!(wZ, Adjoint([13im 14im]), 3, 0.0) &&
                    wZ == setindex!(copy(acomplexintmatmat), hcat([-13im, -14im]), 3, 0.0))
            @test (wZ = QuasiTranspose(copy(complexintmatmat));
                    wZ === setindex!(wZ, Transpose([13im 14im]), 3, 0.0) &&
                    wZ == setindex!(copy(tcomplexintmatmat), hcat([13im, 14im]), 3, 0.0))
        end
    end

    @testset "QuasiAdjoint and Transpose convert methods that convert underlying storage" begin
        intvec, intmat = QuasiVector([1, 2],0:0.5:0.5), QuasiMatrix([1 2 3; 4 5 6],(0:0.5:0.5, Base.OneTo(3)))
        @test convert(QuasiAdjoint{Float64,QuasiVector{Float64,typeof(intvec.axes)}}, QuasiAdjoint(intvec))::QuasiAdjoint{Float64,QuasiVector{Float64,typeof(intvec.axes)}} == QuasiAdjoint(intvec)
        @test convert(QuasiAdjoint{Float64,QuasiMatrix{Float64,typeof(intmat.axes)}}, QuasiAdjoint(intmat))::QuasiAdjoint{Float64,QuasiMatrix{Float64,typeof(intmat.axes)}} == QuasiAdjoint(intmat)
        @test convert(QuasiTranspose{Float64,QuasiVector{Float64,typeof(intvec.axes)}}, QuasiTranspose(intvec))::QuasiTranspose{Float64,QuasiVector{Float64,typeof(intvec.axes)}} == QuasiTranspose(intvec)
        @test convert(QuasiTranspose{Float64,QuasiMatrix{Float64,typeof(intmat.axes)}}, QuasiTranspose(intmat))::QuasiTranspose{Float64,QuasiMatrix{Float64,typeof(intmat.axes)}} == QuasiTranspose(intmat)
    end

    @testset "Adjoint and Transpose similar methods" begin
    intvec, intmat = QuasiVector([1, 2],0:0.5:0.5), QuasiMatrix([1 2 3; 4 5 6],(0:0.5:0.5, Base.OneTo(3)))
        # similar with no additional specifications, vector (rewrapping) semantics
        @test size(similar(QuasiAdjoint(intvec))::QuasiAdjoint{Int,typeof(intvec)}) == size(QuasiAdjoint(intvec))
        @test size(similar(QuasiTranspose(intvec))::QuasiTranspose{Int,typeof(intvec)}) == size(QuasiTranspose(intvec))
        # similar with no additional specifications, matrix (no-rewrapping) semantics
        @test size(similar(QuasiAdjoint(intmat))::QuasiMatrix{Int,typeof(reverse(intmat.axes))}) == size(QuasiAdjoint(intmat))
        @test size(similar(QuasiTranspose(intmat))::QuasiMatrix{Int,typeof(reverse(intmat.axes))}) == size(QuasiTranspose(intmat))
        # similar with element type specification, vector (rewrapping) semantics
        @test size(similar(QuasiAdjoint(intvec), Float64)::QuasiAdjoint{Float64,QuasiVector{Float64,typeof(intvec.axes)}}) == size(QuasiAdjoint(intvec))
        @test size(similar(QuasiTranspose(intvec), Float64)::QuasiTranspose{Float64,QuasiVector{Float64,typeof(intvec.axes)}}) == size(QuasiTranspose(intvec))
        # similar with element type specification, Quasimatrix (no-rewrapping) semantics
        @test size(similar(QuasiAdjoint(intmat), Float64)::QuasiMatrix{Float64,typeof(reverse(intmat.axes))}) == size(QuasiAdjoint(intmat))
        @test size(similar(QuasiTranspose(intmat), Float64)::QuasiMatrix{Float64,typeof(reverse(intmat.axes))}) == size(QuasiTranspose(intmat))
    end

    @testset "Adjoint and Transpose parent methods" begin
        intvec, intmat = QuasiVector([1, 2],0:0.5:0.5), QuasiMatrix([1 2 3; 4 5 6],(0:0.5:0.5, Base.OneTo(3)))
        @test parent(QuasiAdjoint(intvec)) === intvec
        @test parent(QuasiAdjoint(intmat)) === intmat
        @test parent(QuasiTranspose(intvec)) === intvec
        @test parent(QuasiTranspose(intmat)) === intmat
    end

    @testset "Layout" begin
        A = QuasiArray(randn(2,2),([0,0.5],Base.OneTo(2)))
        @test MemoryLayout(typeof(A)) == MemoryLayout(typeof(A')) == MemoryLayout(typeof(transpose(A)))
        @test QuasiArray(A') == QuasiArray(A.parent',(Base.OneTo(2),[0,0.5]))
    end
end


text:
using Mocking

describe_stream_1_patch = @patch function describe_stream(a...; b...)
    return Dict{String, Any}("StreamDescription" => Dict{String, Any}("CreationRequestDateTime" => "Wed May 20 13:51:10 PDT 2015", "TableName" => "Forum", "StreamLabel" => "2015-05-20T20:51:10.252", "StreamViewType" => "NEW_AND_OLD_IMAGES", "Shards" => Any[Dict{String, Any}("ShardId" => "shardId-00000001414562045508-2bac9cd2", "SequenceNumberRange" => Dict{String, Any}("StartingSequenceNumber" => "20500000000000000910398", "EndingSequenceNumber" => "20500000000000000910398")), Dict{String, Any}("ShardId" => "shardId-00000001414576573621-f55eea83", "ParentShardId" => "shardId-00000001414562045508-2bac9cd2", "SequenceNumberRange" => Dict{String, Any}("StartingSequenceNumber" => "820400000000000001192334", "EndingSequenceNumber" => "820400000000000001192334")), Dict{String, Any}("ShardId" => "shardId-00000001414592258131-674fd923", "ParentShardId" => "shardId-00000001414576573621-f55eea83", "SequenceNumberRange" => Dict{String, Any}("StartingSequenceNumber" => "1683700000000000001135967", "EndingSequenceNumber" => "1683700000000000001135967")), Dict{String, Any}("ShardId" => "shardId-00000001414608446368-3a1afbaf", "ParentShardId" => "shardId-00000001414592258131-674fd923", "SequenceNumberRange" => Dict{String, Any}("StartingSequenceNumber" => "2574600000000000000935255"))], "StreamArn" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-20T20:51:10.252", "KeySchema" => Any[Dict{String, Any}("AttributeName" => "ForumName", "KeyType" => "HASH"), Dict{String, Any}("AttributeName" => "Subject", "KeyType" => "RANGE")], "StreamStatus" => "ENABLED"))
end

get_records_1_patch = @patch function get_records(a...; b...)
    return Dict{String, Any}("NextShardIterator" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-20T20:51:10.252|1|AAAAAAAAAAGQBYshYDEe ... <remaining characters omitted> ...", "Records" => Any[Dict{String, Any}("eventVersion" => "1.0", "eventID" => "e2fd9c34eff2d779b297b26f5fef4206", "eventName" => "INSERT", "awsRegion" => "us-west-2", "eventSource" => "aws:dynamodb", "dynamodb" => Dict{String, Any}("SizeBytes" => 41, "SequenceNumber" => "300000000000000499659", "StreamViewType" => "KEYS_ONLY", "Keys" => Dict{String, Any}("ForumName" => Dict{String, Any}("S" => "DynamoDB"), "Subject" => Dict{String, Any}("S" => "DynamoDB Thread 3")), "ApproximateCreationDateTime" => "1.46480646E9")), Dict{String, Any}("eventVersion" => "1.0", "eventID" => "4b25bd0da9a181a155114127e4837252", "eventName" => "MODIFY", "awsRegion" => "us-west-2", "eventSource" => "aws:dynamodb", "dynamodb" => Dict{String, Any}("SizeBytes" => 41, "SequenceNumber" => "400000000000000499660", "StreamViewType" => "KEYS_ONLY", "Keys" => Dict{String, Any}("ForumName" => Dict{String, Any}("S" => "DynamoDB"), "Subject" => Dict{String, Any}("S" => "DynamoDB Thread 1")), "ApproximateCreationDateTime" => "1.46480527E9")), Dict{String, Any}("eventVersion" => "1.0", "eventID" => "740280c73a3df7842edab3548a1b08ad", "eventName" => "REMOVE", "awsRegion" => "us-west-2", "eventSource" => "aws:dynamodb", "dynamodb" => Dict{String, Any}("SizeBytes" => 41, "SequenceNumber" => "500000000000000499661", "StreamViewType" => "KEYS_ONLY", "Keys" => Dict{String, Any}("ForumName" => Dict{String, Any}("S" => "DynamoDB"), "Subject" => Dict{String, Any}("S" => "DynamoDB Thread 2")), "ApproximateCreationDateTime" => "1.46480646E9"))])
end

get_shard_iterator_1_patch = @patch function get_shard_iterator(a...; b...)
    return Dict{String, Any}("ShardIterator" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-20T20:51:10.252|1|AAAAAAAAAAEvJp6D+zaQ...  <remaining characters omitted> ...")
end

list_streams_1_patch = @patch function list_streams(a...; b...)
    return Dict{String, Any}("Streams" => Any[Dict{String, Any}("TableName" => "Forum", "StreamLabel" => "2015-05-20T20:51:10.252", "StreamArn" => "arn:aws:dynamodb:us-wesst-2:111122223333:table/Forum/stream/2015-05-20T20:51:10.252"), Dict{String, Any}("TableName" => "Forum", "StreamLabel" => "2015-05-20T20:50:02.714", "StreamArn" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-20T20:50:02.714"), Dict{String, Any}("TableName" => "Forum", "StreamLabel" => "2015-05-19T23:03:50.641", "StreamArn" => "arn:aws:dynamodb:us-west-2:111122223333:table/Forum/stream/2015-05-19T23:03:50.641")])
end


text:
using LinearAlgebra
using LinearAlgebra.BLAS: BlasFloat, BlasReal   # Union{Complex{Float32}, Complex{Float64}, Float32, Float64}
using LinearAlgebra: checksquare, AbstractTriangular

export BlasFloat, BlasReal,
       BlasVector,BlasMatrix,BlasArray,
       BlasTranspose, BlasAdjoint,
       BlasTransRow, BlasAdjRow, BlasRow,
       BlasNode,
       TransposeTriangular, AdjointTriangular, BlasTriangular, TransformedTriangular,
       uplo_eff, uplochar

BlasVector{T} = AbstractArray{T,1} where  T <: BlasFloat
BlasMatrix{T} = AbstractArray{T,2} where  T <: BlasFloat
BlasArray{T} = Union{BlasVector{T}, BlasMatrix{T}}

BlasTranspose{T} = Transpose{T,P} where P <: BlasArray{T}
BlasAdjoint{T} = Adjoint{T,P} where P <: BlasArray{T}
BlasTransRow{T} = Transpose{T,P} where P <: BlasVector{T}
BlasAdjRow{T} = Adjoint{T,P} where P <: BlasVector{T}
BlasRow{T} = Union{BlasTransRow{T}, BlasAdjRow{T}}

BlasNode{T,N} = Union{Array{T,N}, SubArray{T,N}} where T <:BlasFloat #add more here if needed

TransposeTriangular{T,P} = Transpose{T,P} where P <: AbstractTriangular{T} where T <: BlasFloat
AdjointTriangular{T,P} = Adjoint{T,P} where P <: AbstractTriangular{T} where T <: BlasFloat
TransformedTriangular{T,P} = Union{TransposeTriangular{T,P}, AdjointTriangular{T,P}}

BlasTriangular{T} = Union{ AbstractTriangular{T}, TransformedTriangular{T} } where T <: BlasFloat

AxpyVec{T} = Union{DenseArray{T},AbstractVector{T}} where T <: BlasFloat

transposechar(::AbstractMatrix) = 'N'
transposechar(::Transpose) = 'T'
transposechar(::Adjoint{T}) where T <: Real = 'T'
transposechar(::Adjoint{T}) where T <: Complex = 'C'

blasnode(M::BlasNode) = M  # it's here
blasnode(M::TransformedTriangular) = blasnode(parent(M))
blasnode(M::BlasArray) = blasnode(parent(M))  # go looking for it deeper down ...

uplochar(T::UnitLowerTriangular) = 'L'
uplochar(T::LowerTriangular) = 'L'
uplochar(T::UnitUpperTriangular) = 'U'
uplochar(T::UpperTriangular) = 'U'
uplochar(T::TransformedTriangular) = uplochar(parent(T))

uplo_eff(T::AbstractTriangular) = uplochar(T)
function uplo_eff(T::TransformedTriangular)
    P = uplochar(parent(T))
    P == 'L' && return 'U'
    return 'L'
end

diagchar(T::UnitLowerTriangular) = 'U'
diagchar(T::LowerTriangular) = 'N'
diagchar(T::UnitUpperTriangular) = 'U'
diagchar(T::UpperTriangular) = 'N'
diagchar(T::TransformedTriangular) = diagchar(parent(T))

getparent(T::AbstractTriangular{F}) where F <: BlasFloat = parent(T)
getparent(T::BlasTriangular{F}) where F <: BlasFloat = getfield(T,:parent)

Base.propertynames(T::BlasTriangular) = (:uplo, :diag, :trans, :parent, :blasnode, fieldnames(typeof(T))...)
Base.getproperty(T::BlasTriangular, p::Symbol) = p==:uplo     ? uplochar(T)       :
                                                 p==:diag     ? diagchar(T)       :
                                                 p==:parent   ? getparent(T)      :
                                                 p==:trans    ? transposechar(T)  :
                                                 p==:blasnode ? blasnode(T)       :
                                                 getfield(T,p)
#


text:
"""
    add!!(A, B) -> A′

`A .+= B` if possible; otherwise return `A .+ B`.

# Examples
```jldoctest
julia> using BangBang: add!!

julia> add!!((1,), (2,))
(3,)

julia> add!!([1], [2])
1-element Vector{Int64}:
 3
```
"""
add!!
add!!(A, B) = materialize!!(A, instantiate(broadcasted(+, A, B)))

"""
    mul!!(C, A, B, [α, β]) -> C′
"""
mul!!
mul!!(C, A, B) = may(mul!, C, A, B)
mul!!(C, A, B, α, β) = may(mul!, C, A, B, α, β)

pure(::typeof(mul!)) = NoBang.mul
_asbb(::typeof(mul!)) = mul!!
possible(::typeof(mul!), C, A, B) =
    implements(push!, C) && _matmuleltype(A, B) <: eltype(C)
possible(::typeof(mul!), C, A, B, α, β) =
    implements(push!, C) && _matmuleltype(C, A, B, α, β) <: eltype(C)

# Estimate `eltype` of `C`.  This is how it's done in LinearAlgebra.jl
# but maybe it's better to use the approach of
# https://github.com/tpapp/AlgebraResultTypes.jl ?
_matprod(x, y) = x * y + x * y
_matprod(c, a, b, α, β) = a * b * α + a * b * α + c * β
_matmuleltype(A, B) = Base.promote_op(_matprod, eltype(A), eltype(B))
_matmuleltype(C, A, B, α, β) =
    Base.promote_op(_matprod, eltype(C), eltype(A), eltype(B), typeof(α), typeof(β))

"""
    lmul!!(A, B) -> B′
"""
lmul!!
lmul!!(A, B) = may(lmul!, A, B)

pure(::typeof(lmul!)) = *
_asbb(::typeof(lmul!)) = lmul!!
possible(::typeof(lmul!), A, B) =
    implements(push!, B) && _matmuleltype(A, B) <: eltype(B)

"""
    rmul!!(A, B) -> A′
"""
rmul!!
rmul!!(A, B) = may(rmul!, A, B)

pure(::typeof(rmul!)) = *
_asbb(::typeof(rmul!)) = rmul!!
possible(::typeof(rmul!), A, B) =
    implements(push!, A) && _matmuleltype(A, B) <: eltype(A)


text:
# ENV["CUDA_VISIBLE_DEVICES"] = ""
using Revise
using ADCME
using ADSeismic
using PyPlot
using DelimitedFiles
matplotlib.use("agg")
o = placeholder(ones(1))
param = ElasticPropagatorParams(NX=150, NY=150, NSTEP=5, DELTAT=1e-4, DELTAX=1.0, DELTAY=1.0,
                                vp_ref = 3300.0, USE_PML_XMIN=true, USE_PML_XMAX=true, USE_PML_YMIN=true, USE_PML_YMAX=true)

source = Ricker(param, 15.0, 100.0, 1e6)
srci = repeat([param.NX ÷ 2], 5)
srcj = repeat([param.NY ÷ 2], 5)
srctype = [0;1;2;3;4] # vx, vy, sigmaxx, sigmayy, sigmaxy
srcv = reshape(source, :, 1)
src = ElasticSource(srci, srcj, srctype, srcv)

vp = 3000. * ones(param.NX+2, param.NY+2)
vs = vp / 1.732
rho = 2800. * ones(param.NX+2, param.NY+2)
vp0 = vp
vp_ = Variable(vp[:]) + o[1]^2
vp = reshape(vp_, (152,152))
λ, μ, ρ = compute_lame_parameters(param.NX, param.NY, vp, vs, rho)
model = ElasticPropagatorSolver(param, src, ρ, λ, μ)

sess = Session(); init(sess)

# loss = sum(model.vx)
loss = sum(model.vx^2) + sum(model.vy^2) + sum(model.sigmaxx^2) + sum(model.sigmaxy^2) + sum(model.sigmayy^2)
grad = gradients(loss, o)

# vx = run(sess, model.vx)
# p = visualize_wavefield(vx, param)
# saveanim(p, "elastic-wavefield.gif")

function test(x)

  l = run(sess, loss, o=>x)
  g = run(sess, grad, o=>x)

  return l, g

end

close("all")
ADCME.test_gradients(test, ones(1))
savefig("gradtest_sigmaxy_gpu.png")


text:
const DatasetHandle = Ptr{Void}
const BoosterHandle = Ptr{Void}

const C_API_DTYPE_FLOAT32 = 0
const C_API_DTYPE_FLOAT64 = 1
const C_API_DTYPE_INT32 = 2
const C_API_DTYPE_INT64 = 3

type Dataset
    handle::DatasetHandle

    function Dataset(handle::DatasetHandle)
        ds = new(handle)
        finalizer(ds, Dataset_finalizer)
        return ds
    end

    function Dataset_finalizer(ds::Dataset)
        if ds.handle != C_NULL
            LGBM_DatasetFree(ds)
            ds.handle = C_NULL
        end
    end
end

type Booster
    handle::BoosterHandle
    datasets::Vector{Dataset}

    function Booster(handle::BoosterHandle, datasets::Vector{Dataset})
        bst = new(handle, datasets)
        finalizer(bst, Booster_finalizer)
        return bst
    end

    function Booster_finalizer(bst::Booster)
        if bst.handle != C_NULL
            LGBM_BoosterFree(bst)
            bst.handle = C_NULL
        end
    end
end

function Booster()
    return Booster(C_NULL, Dataset[])
end

function Booster(handle::BoosterHandle)
    return Booster(handle, Dataset[])
end

function jltype_to_lgbmid(datatype::Type)
    if datatype == Float32
        return C_API_DTYPE_FLOAT32
    elseif datatype == Float64
        return C_API_DTYPE_FLOAT64
    elseif datatype == Int32
        return C_API_DTYPE_INT32
    elseif datatype == Int64
        return C_API_DTYPE_INT64
    else
        error("unsupported datatype, got ", datatype)
    end
end

function lgbmid_to_jltype(id::Integer)
    if id == C_API_DTYPE_FLOAT32
        return Float32
    elseif id == C_API_DTYPE_FLOAT64
        return Float64
    elseif id == C_API_DTYPE_INT32
        return Int32
    elseif id == C_API_DTYPE_INT64
        return Int64
    else
        error("unknown LightGBM return type id, got ", id)
    end
end

macro lightgbm(f, params...)
    return quote
        err = ccall(($f, LGBM_library), Cint,
                    ($((esc(i.args[end]) for i in params)...),),
                    $((esc(i.args[end - 1]) for i in params)...))
        if err != 0
            msg = unsafe_string(ccall((:LGBM_GetLastError, LGBM_library), Cstring, ()))
            error("call to LightGBM's ", string($(esc(f))), " failed: ", msg)
        end
    end
end

# function LGBM_DatasetCreateFromFile()
# function LGBM_DatasetCreateFromCSR()
# function LGBM_DatasetCreateFromCSC()

function LGBM_DatasetCreateFromMat{T<:Union{Float32,Float64}}(data::Matrix{T}, parameters::String,
                                                              is_row_major::Bool = false)
    lgbm_data_type = jltype_to_lgbmid(T)
    nrow, ncol = ifelse(is_row_major, reverse(size(data)), size(data))
    out = Ref{DatasetHandle}()
    @lightgbm(:LGBM_DatasetCreateFromMat,
              data => Ptr{Void},
              lgbm_data_type => Cint,
              nrow => Int32,
              ncol => Int32,
              is_row_major => Cint,
              parameters => Cstring,
              C_NULL => Ptr{Void},
              out => Ref{DatasetHandle})
    return Dataset(out[])
end

function LGBM_DatasetCreateFromMat{T<:Real}(data::Matrix{T}, parameters::String)
    return LGBM_DatasetCreateFromMat(convert(Matrix{Float64}, data), parameters)
end

function LGBM_DatasetCreateFromMat{T<:Union{Float32,Float64}}(data::Matrix{T}, parameters::String,
                                                              reference::Dataset,
                                                              is_row_major::Bool = false)
    lgbm_data_type = jltype_to_lgbmid(T)
    nrow, ncol = ifelse(is_row_major, reverse(size(data)), size(data))
    out = Ref{DatasetHandle}()
    @lightgbm(:LGBM_DatasetCreateFromMat,
              data => Ptr{Void},
              lgbm_data_type => Cint,
              nrow => Int32,
              ncol => Int32,
              is_row_major => Cint,
              parameters => Cstring,
              reference.handle => DatasetHandle,
              out => Ref{DatasetHandle})
    return Dataset(out[])
end

function LGBM_DatasetCreateFromMat{T<:Real}(data::Matrix{T}, parameters::String,
                                            reference::Dataset, is_row_major::Bool = false)
    return LGBM_DatasetCreateFromMat(convert(Matrix{Float64}, data), parameters, reference,
                                     is_row_major)
end

# Automatically converts to C's zero-based indices.
function LGBM_DatasetGetSubset(ds::Dataset, used_row_indices::Vector{Int32}, parameters::String)
    num_used_row_indices = length(used_row_indices)
    for idx in 1:num_used_row_indices
        @inbounds used_row_indices[idx] -= 1
    end
    out = Ref{DatasetHandle}()
    @lightgbm(:LGBM_DatasetGetSubset,
              ds.handle => DatasetHandle,
              used_row_indices => Ref{Int32},
              num_used_row_indices => Int32,
              parameters => Cstring,
              out => Ref{DatasetHandle})
    for idx in 1:num_used_row_indices
        @inbounds used_row_indices[idx] += 1
    end
    return Dataset(out[])
end

function LGBM_DatasetGetSubset(ds::Dataset, used_row_indices::Vector{Int64}, parameters::String)
    LGBM_DatasetGetSubset(ds, convert(Vector{Int32}, used_row_indices), parameters)
end

function LGBM_DatasetSetFeatureNames(ds::Dataset, feature_names::Vector{String})
    num_feature_names = length(feature_names)
    @lightgbm(:LGBM_DatasetSetFeatureNames,
              ds.handle => DatasetHandle,
              feature_names => Ref{Cstring},
              num_feature_names => Cint)
end

function LGBM_DatasetGetFeatureNames(ds::Dataset)
    n_features = LGBM_DatasetGetNumFeature(ds)
    feature_names = [Vector{UInt8}(256) for i in 1:n_features]
    num_feature_names = Ref{Cint}()
    @lightgbm(:LGBM_DatasetGetFeatureNames,
              ds.handle => DatasetHandle,
              feature_names => Ref{Ptr{UInt8}},
              num_feature_names => Ref{Cint})

    return [unsafe_string(pointer(feature_name)) for feature_name in feature_names[1:num_feature_names[]]]
end

function LGBM_DatasetFree(ds::Dataset)
    @lightgbm(:LGBM_DatasetFree,
              ds.handle => DatasetHandle)
    return nothing
end

function LGBM_DatasetSaveBinary(ds::Dataset, filename::String)
    @lightgbm(:LGBM_DatasetSaveBinary,
              ds.handle => DatasetHandle,
              filename => Cstring)
    return nothing
end

function _LGBM_DatasetSetField{T<:Union{Float32,Float64,Int32}}(ds::Dataset, field_name::String,
                                                        field_data::Vector{T})
    data_type = jltype_to_lgbmid(T)
    num_element = length(field_data)
    @lightgbm(:LGBM_DatasetSetField,
              ds.handle => DatasetHandle,
              field_name => Cstring,
              field_data => Ptr{Void},
              num_element => Cint,
              data_type => Cint)
    return nothing
end

function LGBM_DatasetSetField(ds::Dataset, field_name::String, field_data::Vector{Float32})
    if field_name == "label" || field_name == "weight"
        _LGBM_DatasetSetField(ds, field_name, field_data)
    else
        _LGBM_DatasetSetField(ds, field_name, convert(Vector{Int32}, field_data))
    end
    return nothing
end

function LGBM_DatasetSetField(ds::Dataset, field_name::String, field_data::Vector{Int32})
    if field_name == "group" || field_name == "group_id"
        _LGBM_DatasetSetField(ds, field_name, field_data)
    else
        _LGBM_DatasetSetField(ds, field_name, convert(Vector{Float32}, field_data))
    end
    return nothing
end

function LGBM_DatasetSetField{T<:Real}(ds::Dataset, field_name::String, field_data::Vector{T})
    if field_name == "label" || field_name == "weight"
        _LGBM_DatasetSetField(ds, field_name, convert(Vector{Float32}, field_data))
    elseif field_name == "init_score"
        _LGBM_DatasetSetField(ds, field_name, convert(Vector{Float64}, field_data))
    elseif field_name == "group" || field_name == "group_id"
        _LGBM_DatasetSetField(ds, field_name, convert(Vector{Int32}, field_data))
    end
    return nothing
end

function LGBM_DatasetGetField(ds::Dataset, field_name::String)
    out_len = Ref{Cint}()
    out_ptr = Ref{Ptr{Void}}()
    out_type = Ref{Cint}()
    @lightgbm(:LGBM_DatasetGetField,
              ds.handle => DatasetHandle,
              field_name => Cstring,
              out_len => Ref{Cint},
              out_ptr => Ref{Ptr{Void}},
              out_type => Ref{Cint})
    jl_out_type = lgbmid_to_jltype(out_type[])
    jl_out_ptr = convert(Ptr{jl_out_type}, out_ptr[])
    return copy(unsafe_wrap(Vector{jl_out_type}, jl_out_ptr, out_len[], false))
end

function LGBM_DatasetGetNumData(ds::Dataset)
    out = Ref{Cint}()
    @lightgbm(:LGBM_DatasetGetNumData,
              ds.handle => DatasetHandle,
              out => Ref{Cint})
    return out[]
end

function LGBM_DatasetGetNumFeature(ds::Dataset)
    out = Ref{Cint}()
    @lightgbm(:LGBM_DatasetGetNumFeature,
              ds.handle => DatasetHandle,
              out => Ref{Cint})
    return out[]
end

function LGBM_BoosterCreate(train_data::Dataset, parameters::String)
    out = Ref{BoosterHandle}()
    @lightgbm(:LGBM_BoosterCreate,
              train_data.handle => DatasetHandle,
              parameters => Cstring,
              out => Ref{BoosterHandle})
    return Booster(out[], [train_data])
end

function LGBM_BoosterCreateFromModelfile(filename::String)
    out_num_iterations = Ref{Cint}()
    out = Ref{BoosterHandle}()
    @lightgbm(:LGBM_BoosterCreateFromModelfile,
              filename => Cstring,
              out_num_iterations => Ref{Cint},
              out => Ref{BoosterHandle})
    return Booster(out[])
end

function LGBM_BoosterLoadModelFromString(model_str::String)
    out_num_iterations = Ref{Cint}()
    out = Ref{BoosterHandle}()
    @lightgbm(:LGBM_BoosterLoadModelFromString,
              model_str => Cstring,
              out_num_iterations => Ref{Cint},
              out => Ref{BoosterHandle})
    return Booster(out[])
end

function LGBM_BoosterFree(bst::Booster)
    @lightgbm(:LGBM_BoosterFree,
              bst.handle => BoosterHandle)
    return nothing
end

function LGBM_BoosterMerge(bst::Booster, other_booster::Booster)
    @lightgbm(:LGBM_BoosterMerge,
              bst.handle => BoosterHandle,
              other_booster.handle => BoosterHandle)
    return nothing
end

function LGBM_BoosterAddValidData(bst::Booster, valid_data::Dataset)
    @lightgbm(:LGBM_BoosterAddValidData,
              bst.handle => BoosterHandle,
              valid_data.handle => DatasetHandle)
    push!(bst.datasets, valid_data)
    return nothing
end

function LGBM_BoosterResetTrainingData(bst::Booster, train_data::Dataset)
    @lightgbm(:LGBM_BoosterResetTrainingData,
              bst.handle => BoosterHandle,
              train_data.handle => DatasetHandle)
    if length(bst.datasets) > 0
        @inbounds bst.datasets[1] = train_data
    else
        bst.datasets = [train_data]
    end
    return nothing
end

function LGBM_BoosterResetParameter(bst::Booster, parameters::String)
    @lightgbm(:LGBM_BoosterResetParameter,
              bst.handle => BoosterHandle,
              parameters => Cstring)
    return nothing
end

function LGBM_BoosterGetNumClasses(bst::Booster)
    out_len = Ref{Cint}()
    @lightgbm(:LGBM_BoosterGetNumClasses,
              bst.handle => BoosterHandle,
              out_len => Ref{Cint})
    return out_len[]
end

function LGBM_BoosterUpdateOneIter(bst::Booster)
    is_finished = Ref{Cint}()
    @lightgbm(:LGBM_BoosterUpdateOneIter,
              bst.handle => BoosterHandle,
              is_finished => Ref{Cint})
    return is_finished[]
end

# function LGBM_BoosterUpdateOneIterCustom()

function LGBM_BoosterRollbackOneIter(bst::Booster)
    @lightgbm(:LGBM_BoosterRollbackOneIter,
              bst.handle => BoosterHandle)
    return nothing
end

function LGBM_BoosterGetCurrentIteration(bst::Booster)
    out_iteration = Ref{Cint}()
    @lightgbm(:LGBM_BoosterGetCurrentIteration,
              bst.handle => BoosterHandle,
              out_iteration => Ref{Cint})
    return out_iteration[]
end

function LGBM_BoosterGetEvalCounts(bst::Booster)
    out_len = Ref{Cint}()
    @lightgbm(:LGBM_BoosterGetEvalCounts,
              bst.handle => BoosterHandle,
              out_len => Ref{Cint})
    return out_len[]
end

function LGBM_BoosterGetEvalNames(bst::Booster)
    out_len = Ref{Cint}()
    n_metrics = LGBM_BoosterGetEvalCounts(bst)
    out_strs = [Vector{UInt8}(256) for i in 1:n_metrics]
    @lightgbm(:LGBM_BoosterGetEvalNames,
              bst.handle => BoosterHandle,
              out_len => Ref{Cint},
              out_strs => Ref{Ptr{UInt8}})
    jl_out_strs = [unsafe_string(pointer(out_str)) for out_str in out_strs[1:out_len[]]]
    return jl_out_strs
end

function LGBM_BoosterGetFeatureNames(bst::Booster)
    out_len = Ref{Cint}()
    n_features = LGBM_BoosterGetNumFeature(bst)
    out_strs = [Vector{UInt8}(256) for i in 1:n_features]
    @lightgbm(:LGBM_BoosterGetFeatureNames,
              bst.handle => BoosterHandle,
              out_len => Ref{Cint},
              out_strs => Ref{Ptr{UInt8}})
    jl_out_strs = [unsafe_string(pointer(out_str)) for out_str in out_strs[1:out_len[]]]
    return jl_out_strs
end

function LGBM_BoosterGetNumFeature(bst::Booster)
    out_len = Ref{Cint}()
    @lightgbm(:LGBM_BoosterGetNumFeature,
              bst.handle => BoosterHandle,
              out_len => Ref{Cint})
    return out_len[]
end

function LGBM_BoosterGetEval(bst::Booster, data::Integer)
    n_metrics = LGBM_BoosterGetEvalCounts(bst)
    out_results = Array{Cdouble}(n_metrics)
    out_len = Ref{Cint}()
    @lightgbm(:LGBM_BoosterGetEval,
              bst.handle => BoosterHandle,
              data => Cint,
              out_len => Ref{Cint},
              out_results => Ref{Cdouble})
    return out_results[1:out_len[]]
end

function LGBM_BoosterGetNumPredict(bst::Booster, data_idx::Integer)
    out_len = Ref{Int64}()
    @lightgbm(:LGBM_BoosterGetNumPredict,
              bst.handle => BoosterHandle,
              data_idx => Cint,
              out_len => Ref{Int64})
    return out_len[]
end

function LGBM_BoosterGetPredict(bst::Booster, data_idx::Integer)
    out_len = Ref{Int64}()
    num_class = LGBM_BoosterGetNumClasses(bst)
    num_data = LGBM_BoosterGetNumPredict(bst, data_idx)
    out_results = Array{Cdouble}(num_class * num_data)
    @lightgbm(:LGBM_BoosterGetPredict,
              bst.handle => BoosterHandle,
              data_idx => Cint,
              out_len => Ref{Int64},
              out_results => Ref{Cdouble})
    return out_results[1:out_len[]]
end

# function LGBM_BoosterPredictForFile()

function LGBM_BoosterCalcNumPredict(bst::Booster, num_row::Integer, predict_type::Integer,
                                    num_iteration::Int)
    out_len = Ref{Int64}()
    @lightgbm(:LGBM_BoosterCalcNumPredict,
              bst.handle => BoosterHandle,
              num_row => Cint,
              predict_type => Cint,
              num_iteration => Cint,
              out_len => Ref{Int64})
    return out_len[]
end

# function LGBM_BoosterPredictForCSR()
# function LGBM_BoosterPredictForCSC()

function LGBM_BoosterPredictForMat{T<:Union{Float32,Float64}}(bst::Booster, data::Matrix{T},
                                                              predict_type::Integer,
                                                              num_iteration::Integer,
                                                              is_row_major::Bool = false)
    num_class = LGBM_BoosterGetNumClasses(bst)
    lgbm_data_type = jltype_to_lgbmid(T)
    nrow, ncol = ifelse(is_row_major, reverse(size(data)), size(data))
    out_len = Ref{Int64}()
    alloc_len = LGBM_BoosterCalcNumPredict(bst, nrow, predict_type, num_iteration)
    out_result = Array{Cdouble}(alloc_len)

    parameter = ""  # full prediction, no early stopping
    @lightgbm(:LGBM_BoosterPredictForMat,
              bst.handle => BoosterHandle,
              data => Ptr{Void},
              lgbm_data_type => Cint,
              nrow => Int32,
              ncol => Int32,
              is_row_major => Cint,
              predict_type => Cint,
              num_iteration => Cint,
              parameter => Cstring,
              out_len => Ref{Int64},
              out_result => Ref{Cdouble})
    return out_result[1:out_len[]]
end

function LGBM_BoosterPredictForMat{T<:Real}(bst::Booster, data::Matrix{T}, predict_type::Integer,
                                            num_iteration::Integer)
    return LGBM_BoosterPredictForMat(bst, convert(Matrix{Float64}, data), predict_type,
                                     num_iteration)
end

function LGBM_BoosterSaveModel(bst::Booster, num_iteration::Integer, filename::String)
    @lightgbm(:LGBM_BoosterSaveModel,
              bst.handle => BoosterHandle,
              num_iteration => Cint,
              filename => Cstring)
    return nothing
end

# function LGBM_BoosterSaveModelToString()
# function LGBM_BoosterDumpModel()
# function LGBM_BoosterGetLeafValue()
# function LGBM_BoosterSetLeafValue()


text:
export MvnModel, mvn_loss_4obj, mvn_loss_1obj

function mvn_logpx(θ::Vector, data::Matrix{Float64})
    n_vars = size(θ, 1) ÷ 2
    μ = θ[1:n_vars]
    Σ = Diagonal(θ[n_vars+1:end] .^ 2)
    dist = MvNormal(μ, Σ)
    logpx = logpdf(dist, data)
    logpx
end

function mvn_loss_4obj(data::Matrix{Float64})
    function p_loss(θ::Vector)
        loss_4obj(mvn_logpx(θ, data))
    end
    p_loss
end

function mvn_loss_1obj(data::Matrix{Float64})
    function p_loss(θ::Vector)
        loss_1obj(mvn_logpx(θ, data))
    end
    p_loss
end

MLJBase.@mlj_model mutable struct MvnModel <: PDFEstimator
    adtype::SciMLBase.AbstractADType = GalacticOptim.AutoZygote()

    optms::Vector{OptM} = long_optms

    loss::Function = mvn_loss_1obj
end

function MLJBase.fit(model::MvnModel, verbosity, X)
    x = collect(MLJBase.matrix(X)')

    n_vars = size(x, 1)
    p = rand(n_vars * 2)
    lss_f = model.loss(x)
    res = optimizeit(model, lss_f, p)

    μ = res.u[1:n_vars]
    Σ = Diagonal(res.u[n_vars+1:end] .^ 2)
    learned_dist = MvNormal(μ, Σ)

    fitresult = (learned_dist, res, n_vars)
    cache = nothing
    report = nothing
    fitresult, cache, report
end

function MLJBase.transform(model::MvnModel, fitresult, Xnew)
    xnew = collect(MLJBase.matrix(Xnew)')
    learned_dist, res, n_vars = fitresult

    ynew = pdf(learned_dist, xnew)
    ynew = reshape(ynew, size(ynew, 1), 1)
    ynew = DataFrame(ynew, :auto)
end

function MLJBase.fitted_params(model::MvnModel, fitresult)
    learned_dist, res, n_vars = fitresult

    (
        learned_params=res.u,
        learned_dist=learned_dist,
        res=res,
        n_vars=n_vars,
    )
end


text:
using GaussianProcesses
using GaussianProcesses: distance, KernelData
import ScikitLearnBase

d, n = 3, 5

x = 2π * rand(d, n)
y = Float64[sum(sin(x[:,i])) for i in 1:n]/d
mZero = MeanZero()
kern = SE(0.0,0.0)

gp = GP(x, y, mZero, kern)

# Function verifies that predictive mean at input observations
# are the same as the output observations
function test_pred_matches_obs(gp::GP)
    y_pred, sig = predict(gp, x)
    @test_approx_eq_eps maximum(abs(gp.y - y_pred)) 0.0 1e-4
end

test_pred_matches_obs(gp)

function sk_test_pred_matches_obs() # ScikitLearn interface test
    gp_sk = ScikitLearnBase.fit!(GP(), x', y)
    y_pred = ScikitLearnBase.predict(gp_sk, x')
    @test_approx_eq_eps maximum(abs(gp_sk.y - y_pred)) 0.0 1e-4
end

sk_test_pred_matches_obs()

# Modify kernel and update
gp.k.ℓ2 = 4.0
x_pred = 2π * rand(d, n)
GaussianProcesses.update_mll!(gp)
y_pred, sig = predict(gp, x_pred)


text:
# docs/make.jl
# Make Documentation for {PROJECT_NAME}.jl

push!(LOAD_PATH, "../src/")

using Documenter, {PROJECT_NAME}

makedocs(
    modules   = [{PROJECT_NAME}],
    doctest   = true,
    clean     = false,
    linkcheck = false,
    format    = Documenter.HTML(prettyurls=!("local" in ARGS)),
    sitename  = "{PROJECT_NAME}.jl",
    authors   = "{FULL_NAME}",
    pages     = [
        "Home" => "index.md",
    ],
)

deploydocs(
    repo = "github.com/{USERNAME}/{PROJECT_NAME}.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)

max_stars_repo_name: JuliaTagBot/julia-template

max_issues_repo_name: JuliaTagBot/julia-template

max_forks_repo_name: JuliaTagBot/julia-template


text:
"""
    Population

Description of the photon population.
"""
mutable struct Population
    n::Int64                           # Size of the population
    r::Vector{SVector{3, Float64}}     # Positions
    μ::Vector{SVector{3, Float64}}     # Directions
    t::Vector{Float64}                 # Time
    isactive::Vector{Bool}             # Is active?
end

"""
    Population(n)

Create an uninitialized population structure with space to hold up to `n`
photons.
"""
Population(n::Int64) = Population(n,
                                  Vector{SVector{3, Float64}}(undef, n),
                                  Vector{SVector{3, Float64}}(undef, n),
                                  Array{Float64, 1}(undef, n),
                                  Array{Bool, 1}(undef, n))

"""
    initphotons!(p, params)

Initialize all photons in a population `p` according to the parameters `params`.
"""
function initphotons!(p::Population, params::Params)
    @unpack source_altitude, source_extension = params

    for i in 1:p.n
        z = source_altitude + (rand() - 0.5) * source_extension
        p.r[i] = @SVector [0, 0, z]
        p.μ[i] = randsphere()
        p.t[i] = 0.0
        p.isactive[i] = true
    end
end


text:
# These benchmarks need to be revised, and run not as part of the tests, but separately.
# https://juliaci.github.io/PkgBenchmark.jl/stable/

#=

"Prints a benchmark results just like btime"
function print_benchmark(trialmin)
    trialtime = BenchmarkTools.time(trialmin)
    trialallocs = BenchmarkTools.allocs(trialmin)
    println(
        "  ",
        BenchmarkTools.prettytime(trialtime),
        " (",
        trialallocs,
        " allocation",
        trialallocs == 1 ? "" : "s",
        ": ",
        BenchmarkTools.prettymemory(BenchmarkTools.memory(trialmin)),
        ")",
    )
end

# test/run_hbv.jl
benchmark = @benchmark Wflow.update(model)
trialmin = BenchmarkTools.minimum(benchmark)
println("HBV Model update")
print_benchmark(trialmin)

# test/run_sbm_gwf.jl
benchmark = @benchmark Wflow.run(tomlpath)
trialmin = BenchmarkTools.minimum(benchmark)
println("SBM GWF Model update (run)")
print_benchmark(trialmin)

# test/run_sbm.jl
benchmark = @benchmark Wflow.update(model)
trialmin = BenchmarkTools.minimum(benchmark)
println("SBM Model update")
print_benchmark(trialmin)

# test/run_sediment.jl
benchmark = @benchmark Wflow.run(tomlpath)
trialmin = BenchmarkTools.minimum(benchmark)
println("Sediment Model update (run)")
print_benchmark(trialmin)
=#


text:
#==============================================================================#
#
#    Declarations of the main `types` used to define Ito diffusions
#    !! IMPORTANT !!
#    to define a diffusion one needs to define a struct and provide a drift
#    and a volatility coefficient functions:
#    `b`: denotes a drift
#    `σ`: denotes a volatility coefficient
#    `a`: (=:σσ') denotes a diffusivity coefficient
#    for linear diffusions
#    `B` and `β` jointly define the drift via `b(t,x):=Bₜx+βₜ`.
#
#==============================================================================#

"""
    DiffusionProcess{T,DP,DW}

Types inheriting from `DiffusionProcess` define Ito diffusions. `T` denotes the
datatype of each coordinate, `DP` the dimension of the stochastic process,
`DW` the dimension of the Wiener process, `SS` lists the state space
restrictions.
"""
abstract type DiffusionProcess{T,DP,DW,SS} end

a(t, x, P::DiffusionProcess) = σ(t, x, P) * σ(t, x, P)'

"""
    dimension(::DiffusionProcess{T,DP,DW})

Return dimension of the stochastic process and driving Brownian motion.
"""
dimension(::DiffusionProcess{T,DP,DW}) where {T,DP,DW} = (
    process = DP,
    wiener = DW
)

"""
    Base.eltype(::DiffusionProcess{T}) where T = T

Return the datatype that each coordinate of the stochastic process is stored in.
"""
Base.eltype(::DiffusionProcess{T}) where T = T

"""
    state_space(::DiffusionProcess{T,DP,DW,SS})

Return the state space restrictions.
"""
state_space(::DiffusionProcess{T,DP,DW,SS}) where {T,DP,DW,SS} = SS

"""
    LinearDiffusion{T,DP,DW,SS} <: DiffusionProcess{T,DP,DW,SS}

Types inheriting from `LinearDiffusion` define a linear Ito-type diffusion, i.e.
solutions to stochastic differential equations of the form:
dXₜ = (BₜXₜ + βₜ)dt + σₜdWₜ, t∈[0,T], X₀=x₀.
"""
abstract type LinearDiffusion{T,DP,DW,SS} <: DiffusionProcess{T,DP,DW,SS} end

"""
    DiffusionDomain

Types inheriting from `DiffusionStateSpace` define the types of restrictions put
on the state-space of the stochastic process.
"""
abstract type DiffusionStateSpace end

const IndexedTime = Tuple{T,S} where {T<:Number,S<:Integer}

#------------------------------------------------------------------------------#
#       Names for functions defining drifts and volatility coefficients
#------------------------------------------------------------------------------#

"""
    b

Compute the drift of a diffusion out-of-place (should use StaticArrays).
"""
function b end

"""
    B

Compute matrix `B` of a linear diffusion out-of-place (should use StaticArrays).
"""
function B end

"""
    β

Compute vector `β` of a linear diffusion out-of-place (should use StaticArrays).
"""
function β end

"""
    σ

Compute the volatility coefficient of a diffusion out-of-place (should use
StaticArrays).
"""
function σ end

"""
    a

Compute the diffusion function σσ' of a diffusion out-of-place (should use
StaticArrays).
"""
function a end

"""
    b!

Compute the drift of a diffusion in-place (uses buffers to store temporary
results).
"""
function b! end

"""
    B!

Compute matrix `B` of a linear diffusion in-place (uses buffers to store
temporary results).
"""
function B! end

"""
    β!

Compute vector `β!` of a linear diffusion in-place (uses buffers to store
temporary results).
"""
function β! end

"""
    σ!

Compute the volatility coefficient of a diffusion in-place (uses buffers to
store temporary results).
"""
function σ! end

"""
    a!

Compute the diffusion function σσ' of a diffusion in-place (uses buffers to
store temporary results).
"""
function a! end

function constdiff end

#------------------------------------------------------------------------------#
#                     Default fallbacks for linear diffusions
#------------------------------------------------------------------------------#
function b!(buffer, t, x, P::LinearDiffusion)
    B!(buffer, t, P)
    β!(buffer, t, P)
    mul!(buffer.b, buffer.B, x, true, true)
end

σ!(buffer, t, x, P::LinearDiffusion) = σ!(buffer, t, P)
a!(buffer, t, x, P::LinearDiffusion) = a!(buffer, t, P)

b(t, x, P::LinearDiffusion) = B(t, P)*x + β(t, P)
σ(t, x, P::LinearDiffusion) = σ(t, P)
a(t, x, P::LinearDiffusion) = a(t, P)
a(t, P::LinearDiffusion) = σ(t, P) * σ(t, P)'

#------------------------------------------------------------------------------#
#
#                   Default fallbacks for indexed time
#
#------------------------------------------------------------------------------#
_b((t,i)::IndexedTime, x, P::DiffusionProcess) = b(t, x, P)
_σ((t,i)::IndexedTime, x, P::DiffusionProcess) = σ(t, x, P)
_b!(buffer, (t,i)::IndexedTime, x, P::DiffusionProcess) = b!(buffer, t, x, P)
_σ!(buffer, (t,i)::IndexedTime, x, P::DiffusionProcess) = σ!(buffer, t, x, P)


text:
import Base: *, transpose, permutedims, A_mul_B!
import Base: A_mul_Bt, A_mul_Bt!, A_mul_Bc, A_mul_Bc!
import Base: At_mul_B, At_mul_B!, Ac_mul_B, Ac_mul_B!
import Base: At_mul_Bt, At_mul_Bt!, Ac_mul_Bc, Ac_mul_Bc!
import Base.LinAlg.BLAS: gemm!
import Base.LinAlg: axpy!, scale!
export axpy!

A_mul_B!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('N','N',one(T),A,B,zero(T),C)
(*){T}(A::KnetMatrix{T},B::KnetMatrix{T})=A_mul_B!(similar(A,(size(A,1),size(B,2))),A,B)

A_mul_Bt!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('N','T',one(T),A,B,zero(T),C)
A_mul_Bt{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt!(similar(A,(size(A,1),size(B,1))),A,B)
A_mul_Bc!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt!(C,A,B)
A_mul_Bc{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=A_mul_Bt(A,B)

At_mul_B!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('T','N',one(T),A,B,zero(T),C)
At_mul_B{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B!(similar(A,(size(A,2),size(B,2))),A,B)
Ac_mul_B!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B!(C,A,B)
Ac_mul_B{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_B(A,B)

At_mul_Bt!{T}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=gemm!('T','T',one(T),A,B,zero(T),C)
At_mul_Bt{T}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt!(similar(A,(size(A,2),size(B,1))),A,B)
Ac_mul_Bc!{T<:Real}(C::KnetMatrix{T}, A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt!(C,A,B)
Ac_mul_Bc{T<:Real}(A::KnetMatrix{T}, B::KnetMatrix{T})=At_mul_Bt(A,B)

function gemm!{T}(transA::Char, transB::Char, alpha::Number, A::KnetArray{T}, B::KnetArray{T}, beta::Number, C::KnetArray{T})
    cublasop(c::Char)=(if c=='N'; 0; elseif c=='T'; 1; elseif c=='C'; 2; else error("Unknown cublas op $c"); end)
    size2(x,i)=(if ndims(x)<=2; size(x,i); elseif i==1; div(length(x),size(x,ndims(x))); elseif i==2; size(x,ndims(x)); else 1; end)
    if transA == 'N'
        m=size2(A,1); k=size2(A,2)
    else
        m=size2(A,2); k=size2(A,1)
    end
    if transB == 'N'
        n=size2(B,2); k==size2(B,1) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    else
        n=size2(B,1); k==size2(B,2) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    end
    (m == size2(C,1) && n == size(C,2)) || throw(DimensionMismatch("$(map(size,(A,B,C)))"))
    transa = cublasop(transA); transb = cublasop(transB)
    alpha = T[alpha]; beta = T[beta]
    lda = size2(A,1); ldb = size2(B,1); ldc = size2(C,1)
    if T<:Float64
        @cuda(cublas, cublasDgemm_v2, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
    elseif T<:Float32
        @cuda(cublas, cublasSgemm_v2, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
    # elseif T<:Float16
    #     @cuda(cublas, cublasHgemm, (Cptr, UInt32, UInt32, Cint, Cint, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle(), transa, transb, m, n, k, alpha, A, lda, B, ldb, beta, C, ldc)
    else
        error("CUBLAS does not support $T")
    end
    return C
end

function axpy!{T}(n::Integer, alpha::Number, x::KnetArray{T}, incx::Integer, y::KnetArray{T}, incy::Integer)
    length(x) == length(y) || throw(DimensionMismatch("$(map(size,(x,y)))"))
    alpha = T[alpha]
    if T<:Float32
        @cuda(cublas, cublasSaxpy_v2, (Cptr, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint), cublashandle(), n, alpha, x, incx, y, incy)
    elseif T<:Float64
        @cuda(cublas, cublasDaxpy_v2, (Cptr, Cint, Ptr{T}, Ptr{T}, Cint, Ptr{T}, Cint), cublashandle(), n, alpha, x, incx, y, incy)
    else
        error("$T not supported")
    end
    return y
end

axpy!{T}(alpha::Number, x::KnetArray{T}, y::KnetArray{T})=axpy!(length(x),alpha,x,1,y,1)

function scal!{T}(n::Integer, alpha::Number, x::KnetArray{T}, incx::Integer)
    alpha = T[alpha]
    if T<:Float32
        @cuda(cublas, cublasSscal_v2, (Cptr, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle(), n, alpha, x, incx)
    elseif T<:Float64
        @cuda(cublas, cublasDscal_v2, (Cptr, Cint, Ptr{T}, Ptr{T}, Cint), cublashandle(), n, alpha, x, incx)
    else
        error("$T not supported")
    end
    return x
end

scale!{T}(alpha::Number, x::KnetArray{T})=scal!(length(x),alpha,x,1)
scale!{T}(x::KnetArray{T}, alpha::Number)=scal!(length(x),alpha,x,1)

function transpose{T}(x::KnetArray{T})
    ndims(x) != 2 && error("Transpose is supported only for 2D KnetArrays")
    sz = size(x)
    y = similar(x,(sz[2],sz[1]))
    if T<:Float32
        @cuda(cublas, cublasSgeam, (Cptr,UInt32,UInt32,Cint,Cint,Ptr{T},Ptr{T},Cint,Ptr{T},Ptr{T},Cint,Ptr{T},Cint),
              cublashandle(),1,1,size(y,1),size(y,2),Ref(T(1.0)),x,size(x,1),Ref(T(0.0)),x,size(x,1),y,size(y,1))
    elseif T<:Float64
        @cuda(cublas, cublasDgeam, (Cptr,UInt32,UInt32,Cint,Cint,Ptr{T},Ptr{T},Cint,Ptr{T},Ptr{T},Cint,Ptr{T},Cint),
              cublashandle(),1,1,size(y,1),size(y,2),Ref(T(1.0)),x,size(x,1),Ref(T(0.0)),x,size(x,1),y,size(y,1))
    else
        error("CUBLAS does not support $T")
    end
    return y
end

"""

    mat(x)

Reshape x into a two-dimensional matrix.

This is typically used when turning the output of a 4-D convolution
result into a 2-D input for a fully connected layer.  For 1-D inputs
returns `reshape(x, (length(x),1))`.  For inputs with more than two
dimensions of size `(X1,X2,...,XD)`, returns

    reshape(x, (X1*X2*...*X[D-1],XD))

"""
function mat(x)
    if ndims(x) > 2
        xn = size(x,ndims(x))
        reshape(x, (div(length(x),xn),xn))
    elseif ndims(x)==2
        x
    elseif ndims(x)==1
        reshape(x, (length(x),1))
    else
        throw(MethodError(mat,x))
    end
end

import Base: permutedims, ipermutedims
function permutedims{T,N}(x::KnetArray{T,N}, dims)
    if length(dims) != N; throw(DimensionMismatch()); end
    if N == 2
        # use individual == to cover row col vecs tuples etc
        if dims[1]==1 && dims[2]==2
            return copy(x)
        elseif dims[1]==2 && dims[2]==1
            return transpose(x)
        else
            throw(ArgumentError("no valid permutation of dimensions"))
        end
    elseif N == 3
        if dims[1]==1 && dims[2]==2 && dims[3]==3
            return copy(x)
        else
            funcName = permutefunc(x,dims)
            y = similar(x, size(x,dims[1]), size(x,dims[2]), size(x,dims[3]))
            @eval ccall(($funcName,libknet8),Void,(Ptr{$T},Cint,Cint,Cint,Ptr{$T},Cint,Cint,Cint),
                        $x,size($x,1),size($x,2),size($x,3),$y,size($y,1),size($y,2),size($y,3))
            return y
        end
    elseif N == 4 || N == 5
        error("Not yet implemented")
    end
end

function permutefunc{T,N}(x::KnetArray{T,N}, dims)
    funcName = "permutedims$(N)D_"
    for i=1:N
        funcName = funcName * "$(dims[i])_"
    end
    if T<:Float32
        funcName = funcName * "32"
    elseif T<:Float64
        funcName = funcName * "64"
    else
        error("$T not supported")
    end
    funcName = funcName * "_44"
    return funcName
end

function ipermutedims(A::KnetArray,perm)
    iperm = Array{Int}(length(perm))
    for (i,p) = enumerate(perm)
        iperm[p] = i
    end
    return permutedims(A,iperm)
end

# Low level gemm! call with pointers

using Base.LinAlg
using Base.LinAlg.BLAS: libblas, BlasInt
using Compat: @blasfunc

# C := alpha*op(A)*op(B) + beta*C, where:
# op(X) is one of op(X) = X, or op(X) = XT, or op(X) = XH,
# alpha and beta are scalars,
# A, B and C are matrices:
# op(A) is an m-by-k matrix,
# op(B) is a k-by-n matrix,
# C is an m-by-n matrix.

for (gemm, elty) in ((:dgemm_,:Float64), (:sgemm_,:Float32))
    @eval begin
        function gemm!(transA::Char, transB::Char, M::Int, N::Int, K::Int, alpha::($elty), A::Ptr{$elty}, B::Ptr{$elty}, beta::($elty), C::Ptr{$elty})
            if transA=='N'; lda=M; else; lda=K; end
            if transB=='N'; ldb=K; else; ldb=N; end
            ldc = M;
            ccall((@blasfunc($gemm), libblas), Void,
                  (Ptr{UInt8}, Ptr{UInt8}, Ptr{BlasInt}, Ptr{BlasInt},
                   Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
                   Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty},
                   Ptr{BlasInt}),
                  &transA, &transB, &M, &N, &K,
                  &alpha, A, &lda, B, &ldb, &beta, C, &ldc)
        end
    end
end