text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* File containg models of columns: distillation, stripping, absorbers
#* rectifier, ....
#*
#* The default nomenclature is:
#*		Type_Column_reboilertype_condensertyper
#*
#* where:
#*	Type = refluxed or reboiled or section
#*	Column = Stripping, Absorption, Rectifier, Distillation
#*	Reboiler type (if exists) = kettle or thermosyphon
#*	Condenser type (if exists) = with subccoling or without subcooling
#*
#*-----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id: column.mso 511 2008-05-12 17:25:33Z paula $
#*---------------------------------------------------------------------
type Packed_Section_ColumnBasic
	Packed_Section_ColumnBasic()=begin
		PP=outers.PP
		NComp=outers.NComp
		[

		],
		[

		],
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of components",
				:Protected=>true
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Feed Stages",
				:Default=>3,
				:Protected=>true
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Stages",
				:Default=>3
			)),
			fill(DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Feed Stage Location",
				:Default=>2
			)),(NumberOfFeeds)),
			length(Dict{Symbol,Any}(
				:Brief=>"Height of packing"
			)),
			heat_rate(Dict{Symbol,Any}(
				:Brief=>"Rate of heat supply"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Column diameter"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Void fraction of packing, (m^3 void space/m^3 packed bed)"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Resistance coefficient on the liquid load",
				:Default=>1
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"surface area per packing volume",
				:Unit=>"m^2/m^3"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"The Height Equivalent to a Theoretical Plate",
				:Protected=>true
			)),
			[
				:(LiquidConnector.F= LiquidInlet.F),
				:(LiquidConnector.T = LiquidInlet.T),
				:(LiquidConnector.P = LiquidInlet.P),
				:(LiquidConnector.z = LiquidInlet.z),
				:(LiquidConnector.v = LiquidInlet.v),
				:(LiquidConnector.h = LiquidInlet.h),
				:(VapourConnector.F= VapourInlet.F),
				:(VapourConnector.T = VapourInlet.T),
				:(VapourConnector.P = VapourInlet.P),
				:(VapourConnector.z = VapourInlet.z),
				:(VapourConnector.v = VapourInlet.v),
				:(VapourConnector.h = VapourInlet.h),
				:(LiquidOutlet.F= STAGES(NumberOfStages).OutletLiquid.F),
				:(LiquidOutlet.T = STAGES(NumberOfStages).OutletLiquid.T),
				:(LiquidOutlet.P = STAGES(NumberOfStages).OutletLiquid.P),
				:(LiquidOutlet.z = STAGES(NumberOfStages).OutletLiquid.z),
				:(VapourOutlet.F= STAGES(1).OutletVapour.F),
				:(VapourOutlet.T = STAGES(1).OutletVapour.T),
				:(VapourOutlet.P = STAGES(1).OutletVapour.P),
				:(VapourOutlet.z = STAGES(1).OutletVapour.z),
				:(TCI*"K" = STAGES(CONTROL.Tindicator_TrayNumber).OutletVapour.T),
				:(PCI*"atm" = STAGES(CONTROL.Pindicator_TrayNumber).OutletVapour.P),
				:(zCI * 1e-6 = STAGES(CONTROL.zindicator_TrayNumber).OutletLiquid.z(CONTROL.zindicator_Component) * STAGES(CONTROL.zindicator_TrayNumber).Mw(CONTROL.zindicator_Component) / sum(STAGES(CONTROL.zindicator_TrayNumber).Mw*STAGES(CONTROL.zindicator_TrayNumber).OutletLiquid.z)),
				:(PressureDrop = STAGES(NumberOfStages).OutletLiquid.P - STAGES(1).OutletLiquid.P),
			],
			[
				"","","","","","","","","","","","","","","","","","","","","Tray Temperature Indicator","Tray Pressure Indicator","Tray Composition Indicator","PressureDrop",
			],
			[:PP,:NComp,:NumberOfFeeds,:NumberOfStages,:FeedStageLocation,:PackingHeight,:HeatSupply,:ColumnDiameter,:VoidFraction,:ResistanceCoeff,:AreaPerPackingVol,:HETP,],
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	NumberOfFeeds::DanaInteger
	NumberOfStages::DanaInteger
	FeedStageLocation::Array{DanaInteger}
	PackingHeight::length
	HeatSupply::heat_rate
	ColumnDiameter::length
	VoidFraction::DanaReal
	ResistanceCoeff::positive
	AreaPerPackingVol::DanaReal
	HETP::length
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	parameters::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Packed_Section_ColumnBasic
function set(in::Packed_Section_ColumnBasic)
	STAGES.Number_Stages = NumberOfStages
	 STAGES.HeightOfPacking = PackingHeight
	 STAGES.HeatOnStage = HeatSupply
	 STAGES.ColumnInternalDiameter = ColumnDiameter
	 STAGES.PackingVoidFraction = VoidFraction
	 STAGES.LiquidResistanceCoeff = ResistanceCoeff
	 STAGES.AreaPerPackingVolume = AreaPerPackingVol
	 HETP = PackingHeight/NumberOfStages

end
function setEquationFlow(in::Packed_Section_ColumnBasic)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	#TCI*'K' = TRAYS(min([NumberOfTrays, CONTROL.Tindicator_TrayNumber])).OutletVapour.T;
	addEquation(21)
	#PCI*'atm' = TRAYS(min([NumberOfTrays, CONTROL.Pindicator_TrayNumber])).OutletVapour.P;
	addEquation(22)
	#zCI * .000001 = STAGES(CONTROL.zindicator_TrayNumber).OutletVapour.z(CONTROL.zindicator_Component) * STAGES(CONTROL.zindicator_TrayNumber).Mw(CONTROL.zindicator_Component) ; #/sum(STAGES(CONTROL.zindicator_TrayNumber).Mw*STAGES(CONTROL.zindicator_TrayNumber).OutletLiquid.z);
	addEquation(23)
	addEquation(24)
end
function initial(in::Packed_Section_ColumnBasic)

end
function atributes(in::Packed_Section_ColumnBasic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/PackedSectionColumn"
	fields[:Brief]="Model of a packed column section."
	fields[:Info]="== Model of a packed column section containing ==
* NStages theoretical stages.

== Specify ==
* the feed stream of each tray (Inlet);
* the InletLiquid stream of the top tray;
* the InletVapour stream of the bottom tray;
* the total pressure drop (dP) of the section.

== Initial Conditions ==
* the stages temperature (OutletLiquid.T);
* the stages liquid holdup;
* (NoComps - 1) OutletLiquid (OR OutletVapour) compositions for each tray.
"
	drive!(fields,_)
	return fields
end
Packed_Section_ColumnBasic(_::Dict{Symbol,Any})=begin
	newModel=Packed_Section_ColumnBasic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
@testset "nl_convert" begin
    @testset "check special conversion cases" begin
        special_cases = [:cbrt, :abs2, :inv, :log2, :log1p, :exp2, :expm1, :sec,
                         :csc, :cot, :sind, :cosd, :tand, :asind, :acosd, :atand,
                         :secd, :cscd, :cotd, :sech, :csch, :coth, :asech, :acsch]
        for func in special_cases
            x = rand()
            expr = Expr(:call, func, x)
            @test isapprox(eval(JAMSDWriter.convert_formula(expr)), eval(expr), atol=1e-6)
        end
        # These functions need input >1
        for func in [:acoth, :asec, :acsc, :acot, :asecd, :acscd, :acotd]
            x = rand() + 1
            expr = Expr(:call, func, x)
            @test isapprox(eval(JAMSDWriter.convert_formula(expr)), eval(expr), atol=1e-6)
        end
    end

    @testset "check numeric values" begin
        x = rand()
        @test JAMSDWriter.convert_formula(:($x)) == :($x)
        x = -rand()
        @test JAMSDWriter.convert_formula(:($x)) == :($x)
    end

    @testset "check unary, binary and n-ary plus" begin
        expr = :(+(1))
        @test JAMSDWriter.convert_formula(expr) == :(1)
        expr = :(1 + 2)
        @test JAMSDWriter.convert_formula(expr) == :(1 + 2)
        expr = :(1 + 2 + 3)
        @test JAMSDWriter.convert_formula(expr) == :(sum(1, 2, 3))
    end

    @testset "check unary, binary and n-ary minus" begin
        expr = :(- x)
        @test JAMSDWriter.convert_formula(expr) == :(neg(x))
        expr = :(x - y)
        @test JAMSDWriter.convert_formula(expr) == :(x - y)
        expr = :(x - y - z)
        @test JAMSDWriter.convert_formula(expr) == :((x - y) - z)
    end

    @testset "check n-ary multiplication" begin
        expr = :(x * y * z)
        @test JAMSDWriter.convert_formula(expr) == :(x * (y * z))
    end

    @testset "check comparison expansion" begin
        expr = :(1 < 2 < 3)
        @test JAMSDWriter.convert_formula(expr) == :(1 < 2 && 2 < 3)
        expr = :(1 < 2 < 3 < 4)
        @test JAMSDWriter.convert_formula(expr) == :((1 < 2 && 2 < 3) && 3 < 4)
    end
end

ext: jl


text:
# Install MySQL
julia> Pkg.add("MySQL")
julia> using MySQL

# Create a database in mySQL
using MySQL
con = mysql_connect("localhost", "username", "password", "db_name")
command = """CREATE TABLE Employee
             (
                 ID INT NOT NULL AUTO_INCREMENT,
                 Name VARCHAR(255),
                 Salary FLOAT,
                 JoinDate DATE,
                 PRIMARY KEY (ID)
             );"""
mysql_execute(con, command)

# Insert some values
mysql_execute(con, "INSERT INTO Employee (Name, Salary, JoinDate) values
('Alpha', 25000.00, '2015-12-12'), ('Beta', 35000.00, '2012-18-17), ('Gamma', 50000.00, '2013-12-14');")

# Get SELECT results
command = "SELECT * FROM Employee;"
dframe = mysql_execute(con, command)

# Close connection
mysql_disconnect(con)

### NoSQL
julia> using Mongo, LibBSON

# Create a client connection
julia> client = MongoClient()

# Get a handle to collection named "test" in database "test".
# Client object, database name, and collection name are stored as variables.
test = MongoCollection(client, "test", "test")

# Insert a document
document = insert(test, Dict("name" => "rahul", "book" => "learning julia"))

### Using JSON
julia> using Requests
julia> using JSON

julia> result = Requests.get("http://httpbin.org/get")
Response(200 OK, 14 headers, 281 bytes in body)

julia> typeof(result)
HttpCommon.Response

julia> fieldnames(result)
8-element Array{Symbol,1}:
 :status
 :headers
 :cookies
 :data
 :request
 :history
 :finished
 :requests

julia> println(result.data)
UInt8[0x7b, 0x0a, 0x20, 0x20, 0x22, 0x61, 0x72, 0x67, 0x73, 0x22, 0x3a, 0x20, 0x7b, 0x7d, 0x2c, 0x20, 0x0a, 0x20, 0x20, 0x22, 0x68, 0x65, 0x61, 0x64, 0x65, 0x72, 0x73, 0x22, 0x3a, 0x20, 0x7b, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x22, 0x41, 0x63, 0x63, 0x65, 0x70, 0x74, 0x22, 0x3a, 0x20, 0x22, 0x74, 0x65, 0x78, 0x74, 0x2f, 0x68, 0x74, 0x6d, 0x6c, 0x2c, 0x61, 0x70, 0x70, 0x6c, 0x69, 0x63, 0x61, 0x74, 0x69, 0x6f, 0x6e, 0x2f, 0x78, 0x68, 0x74, 0x6d, 0x6c, 0x2b, 0x78, 0x6d, 0x6c, 0x2c, 0x61, 0x70, 0x70, 0x6c, 0x69, 0x63, 0x61, 0x74, 0x69, 0x6f, 0x6e, 0x2f, 0x78, 0x6d, 0x6c, 0x3b, 0x71, 0x3d, 0x30, 0x2e, 0x39, 0x2c, 0x2a, 0x2f, 0x2a, 0x3b, 0x71, 0x3d, 0x30, 0x2e, 0x38, 0x22, 0x2c, 0x20, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x22, 0x43, 0x6f, 0x6e, 0x6e, 0x65, 0x63, 0x74, 0x69, 0x6f, 0x6e, 0x22, 0x3a, 0x20, 0x22, 0x63, 0x6c, 0x6f, 0x73, 0x65, 0x22, 0x2c, 0x20, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x22, 0x48, 0x6f, 0x73, 0x74, 0x22, 0x3a, 0x20, 0x22, 0x68, 0x74, 0x74, 0x70, 0x62, 0x69, 0x6e, 0x2e, 0x6f, 0x72, 0x67, 0x22, 0x2c, 0x20, 0x0a, 0x20, 0x20, 0x20, 0x20, 0x22, 0x55, 0x73, 0x65, 0x72, 0x2d, 0x41, 0x67, 0x65, 0x6e, 0x74, 0x22, 0x3a, 0x20, 0x22, 0x52, 0x65, 0x71, 0x75, 0x65, 0x73, 0x74, 0x73, 0x2e, 0x6a, 0x6c, 0x2f, 0x30, 0x2e, 0x30, 0x2e, 0x30, 0x22, 0x0a, 0x20, 0x20, 0x7d, 0x2c, 0x20, 0x0a, 0x20, 0x20, 0x22, 0x6f, 0x72, 0x69, 0x67, 0x69, 0x6e, 0x22, 0x3a, 0x20, 0x22, 0x33, 0x35, 0x2e, 0x31, 0x38, 0x38, 0x2e, 0x31, 0x34, 0x32, 0x2e, 0x31, 0x36, 0x32, 0x22, 0x2c, 0x20, 0x0a, 0x20, 0x20, 0x22, 0x75, 0x72, 0x6c, 0x22, 0x3a, 0x20, 0x22, 0x68, 0x74, 0x74, 0x70, 0x3a, 0x2f, 0x2f, 0x68, 0x74, 0x74, 0x70, 0x62, 0x69, 0x6e, 0x2e, 0x6f, 0x72, 0x67, 0x2f, 0x67, 0x65, 0x74, 0x22, 0x0a, 0x7d, 0x0a]

julia> typeof(result.data)
Array{UInt8,1}

julia> output = JSON.parse(convert(String, result.data))
Dict{String,Any} with 4 entries:
  "headers" => Dict{String,Any}(Pair{String,Any}("Connection", "close"),Pair{St…
  "args"    => Dict{String,Any}()
  "url"     => "http://httpbin.org/get"
  "origin"  => "35.188.142.162"

julia> output["headers"]
Dict{String,Any} with 4 entries:
  "Connection" => "close"
  "Host"       => "httpbin.org"
  "Accept"     => "text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=…
  "User-Agent" => "Requests.jl/0.0.0"

julia> result = Requests.post("http://httpbin.org/post"; data = "this is julia")
Response(200 OK, 14 headers, 436 bytes in body)

julia> JSON.parse(convert(String, result.data))
Dict{String,Any} with 8 entries:
  "headers" => Dict{String,Any}(Pair{String,Any}("Connection", "close"),Pair{St…
  "json"    => nothing
  "files"   => Dict{String,Any}()
  "args"    => Dict{String,Any}()
  "data"    => "this is julia"
  "url"     => "http://httpbin.org/post"
  "form"    => Dict{String,Any}()
  "origin"  => "35.188.142.162"

max_stars_repo_name: PacktPublishing/Learning-Julia

max_issues_repo_name: PacktPublishing/Learning-Julia

max_forks_repo_name: PacktPublishing/Learning-Julia


text:
# type for initial manpower described as [Academic Level, Seniority, Personnel Number]
type InitMPcluster
    Academiclvl::String
    Affiliation::String
    Seniority::Int
    Nb::Int
    ActualSeniority::Int
    InitMPcluster()=new("", 0, 0)
    InitMPcluster(acdlvl::String, affil::String, senior::Int, nb::Int) = new(acdlvl, affil, senior, nb, senior)
end

type ManpowerObjective
    priority::Int
    targetLevel::Vector{DataType}
    Objectives::Vector{String}
    Number::Int
    Relation::String #relation between the diferent levels and / or
    InitTolerance::Float64
    EndTolerance::Float64
    Alfa::Float64
    PlotNb::Int
    ManpowerObjective(pri::Int;
                      target::Vector{DataType} = Vector{DataType}(),
                      obj::Vector{String} = Vector{String}(),
                      nb::Int = 0,
                      rl = "",
                      initT = 0.0,
                      endT = 0.0,
                      alfa = 1.0,
                      plot = 0) = new(pri, target, obj, nb, rl, initT, endT, alfa, plot)
end

TypesDict = Dict{String, DataType}("Academ"=> AcademicLevel,
                                   "Affil"=> Affiliation,
                                   "Assign"=> Assignment,
                                   "Rank"=> Rank,
                                   "Job"=> Job,
                                   "SubJob"=> SubJob)

type Subpopulation
    priority::Int
    Name::String
    NbRequired::Int

    Subpopulation(;pr = 1, nm = "", nb = 0) = new(pr, nm, nb)

end

type Dependency

    Sets::Vector{Vector{Int}}
    Percentages::Vector{Float32}
    Priority::Int

    Dependency(;pr = 0) = new(Vector{Vector{Int}}(), Vector{Float32}(), pr)
end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
if !isdefined(Main, :GenerateData)
  function GenerateData(kappa,omega,Nsubj,Ntrials)
    data = fill(0,Nsubj*Ntrials)
    SubjIdx = similar(data)
    cnt = 0
    for subj in 1:Nsubj
      alpha = kappa*omega
      beta = (1-kappa)*omega
      theta = rand(Beta(alpha,beta))
      for trial in 1:Ntrials
        cnt += 1
        data[cnt] = rand(Bernoulli(theta))
        SubjIdx[cnt] = subj
      end
    end
    return data,SubjIdx
  end
end

if !isdefined(Main, :InitializeChain)
  function InitializeChain(Nteams,Nchains,teamParm,divParm,teamVal,divVal)
    init = Dict{String, Any}[Dict() for chains = 1:Nchains]
    for chain in 1:Nchains
      for team in 1:Nteams
        for i in zip(teamParm,teamVal)
           name = string(i[1],"[",team,"]")
           init[chain][name] = rand(Normal(i[2][1],i[2][2]))
        end
      end
      for i in zip(divParm,divVal)
        name = string(i[1])
        init[chain][name] = i[2]
      end
    end
    return init
  end
end


text:
Base.@kwdef mutable struct GridStruct
    nz = nothing
    zP = nothing
    z = nothing
    dz = nothing
    NumFaces = nothing
    Faces = nothing
    NumEdges = nothing
    Edges = nothing
    NumNodes = nothing
    Nodes = nothing
    Form = nothing
    Type = nothing
    nBar = nothing
    Dim = nothing
    Rad = nothing
    NumEdgesI = nothing
    NumEdgesB = nothing
    vtkP = nothing
    ConnectivityList = nothing
    nBar3 = nothing
    KiteGrid = nothing
end


text:
#include("ConeMantle.jl")
#include("EllipsoidMantle.jl")
#include("EllipticalSurface.jl")
#include("Plane.jl")
#include("Polygon.jl")
#include("TorusMantle.jl")

@recipe function f(s::AbstractSurfacePrimitive{T}; n_arc = 40, n_vert_lines = 2, n_samples = 40) where {T}
    seriestype --> :csg
    l = get_label_name(s)
    if haskey(plotattributes, :seriestype)
        if plotattributes[:seriestype] == :csg
            @series begin
                label := ""
                linewidth := 0
                linecolor --> :white
                mesh(s, n_arc)
            end
            @series begin
                label --> l
                linecolor --> :black
                fillalpha := 1
                linewidth --> 2
                mesh(s, n_arc, n_vert_lines)
            end
        elseif plotattributes[:seriestype] == :mesh3d
            label --> l
            linecolor --> :white
            mesh(s, n_arc)
        elseif plotattributes[:seriestype] == :wireframe
            label --> l
            seriescolor --> :black
            fillalpha := 1
            linewidth --> 2
            mesh(s, n_arc, n_vert_lines)
        elseif plotattributes[:seriestype] == :samplesurface
            label --> l
            seriesalpha --> 0.4
            sample(s, extremum(s)/n_samples)
        else
            @warn "The only seriestypes which will return a plot are :csg, :wireframe, :mesh3d, and :samplesurface"
        end
    end
end

@recipe function f(::Type{Val{:samplesurface}}, x, y, z)
    seriescolor --> 1
    seriesalpha --> 0.2
    markerstrokewidth --> 0
    markersize --> 4
    seriestype := :scatter
    ()
end

@recipe function f(vp::AbstractVector{<:AbstractSurfacePrimitive})
    linecolor --> :black
    @series begin
        label --> "Faces"
        show_normal --> false
        vp[1]
    end
    if length(vp) > 1
        for p in vp[2:end]
            @series begin
                show_normal --> false
                label := nothing
                p
            end
        end
    end
end


text:
function lll_basis_profile(A::NfOrdIdl; prec::Int = 100)
  c = Hecke.minkowski_mat(Hecke.nf(order(A)), prec)
  l = lll(basis_mat(A))
  b = FakeFmpqMat(l)*basis_mat(order(A))
  if !isdefined(rt_c, :cache)
    rt_c.cache = 0*c
  end
  d = rt_c.cache
  Hecke.mult!(d, b.num, c)

  old = get_bigfloat_precision()
  set_bigfloat_precision(prec)

  g = Hecke.round_scale(d, prec)
  set_bigfloat_precision(old)
  g = g*g'
  Hecke.shift!(g, -prec)
  g += rows(g)*one(parent(g))

  l = lll_gram(g)

  lp = [ div(l[i,i], fmpz(2)^prec) for i=1:rows(l)]
  return lp
end

#function short_elem(c::roots_ctx, A::NfOrdIdl,
#                v::fmpz_mat = MatrixSpace(FlintZZ, 1,1)(); prec::Int = 100)
#  l, t = lll(c, A, v, prec = prec)
#  w = window(t, 1,1, 1, cols(t))
#  c = w*b
#  q = elem_from_mat_row(K, c, 1, b_den)
#  return q
#end

function bkz_basis(A::NfOrdIdl, bs::Int;
                      v::fmpz_mat = zero_matrix(FlintZZ, 1, 1),
                      prec::Int = 100)

  K = nf(order(A))

  c = minkowski_mat(K, prec)

  l, t1 = lll_with_transform(basis_mat(A))
  temp = FakeFmpqMat(l)*basis_mat(order(A))
  b = temp.num
  b_den = temp.den

  if !isdefined(rt_c, :cache)
    rt_c.cache = 0*c
  end
  d = rt_c.cache
  mult!(d, b, c)

  #ignore v

  g = round_scale(d, prec)

  #println("calling bkz")
  g, tt = bkz_trans(g, bs)  ## check: bkz_trans seems to s.t. kill the input

  c = tt*b
  q = nf_elem[elem_from_mat_row(K, c, i, b_den) for i=1:degree(K)]

  return q
end

function fplll_basis(rt_c::Hecke.roots_ctx, A::NfOrdIdl, bs::Int;
                      v::fmpz_mat = zero_matrix(FlintZZ, 1,1),
                      prec::Int = 100)

  K = nf(order(A))

  c = minkowski_mat(rt_c, K, prec)

  @time l, t1 = lll_with_transform(basis_mat(A))
  @time temp = FakeFmpqMat(l)*basis_mat(order(A))
  b = temp.num
  b_den = temp.den

  if !isdefined(rt_c, :cache)
    rt_c.cache = 0*c
  end
  d = rt_c.cache
  @time mult!(d, b, c)

  #ignore v

  @time g = round_scale(d, prec)

  println("calling lll")
  @time g, tt = FPlll.lll_trans(g)  ## check: bkz_trans seems to s.t. kill the input

  c = tt*b
  q = nf_elem[elem_from_mat_row(K, c, i, b_den) for i=1:degree(K)]

  return q
end

function random_ideal_with_norm_up_to(a::Hecke.NfFactorBase, B::Integer)
  B = fmpz(B)
  O = order(a.ideals[1])
  K = Hecke.nf(O)
  I = Hecke.ideal(O, K(1))
  while B >= norm(a.ideals[end])
    J = a.ideals[rand(find(x -> (norm(x) <= B), a.ideals))]
    B = div(B, norm(J))
    I = I*J
  end
  return I
end

##chebychev_u ... in flint
function tschebyschew(Qx::Nemo.FmpqPolyRing, n::Int)
  T = [Qx(1), gen(Qx)]
  while length(T) <= n
    push!(T, 2*T[2]*T[end] - T[end-1])
  end
  return T[end]
end

function auto_of_maximal_real(K::AnticNumberField, n::Int)
  # image of zeta -> zeta^n
  # assumes K = Q(zeta+1/zeta)
  # T = tschebyschew(n), then
  # cos(nx) = T(cos(x))
  # zeta + 1/zeta = 2 cos(2pi/n)
  T = tschebyschew(parent(K.pol), n)
  i = evaluate(T, gen(K)*1//fmpz(2))*2
  return Mor(K, K, i)
end

function auto_simplify(A::Map, K::AnticNumberField)
  Qx = parent(K.pol)
  b = A(gen(K))
  return Mor(K, K, b)
end

function auto_power(A::Map, n::Int)
  if n==1
    return A
  end;
  B = x->A(A(x));
  C = auto_power(B, div(n, 2))
  if n%2==0
    return C
  else
    return x-> A(C(x))
  end
end

function orbit(f::Map, a::Nemo.nf_elem)
  b = Set([a])
  nb = 1
  while true
    b = union(Set([f(x) for x in b]) , b)
    if nb == length(b)
      return b
    end
    nb = length(b)
  end
end

function order_of_auto(f::Map, K::AnticNumberField)
  a = gen(K)
  b = f(a)
  i = 1
  while b != a
    b = f(b)
    i += 1
  end
  return i
end


text:
## -------------------------------------------------------------------
function load_model(modelkey::String; uncompress = true)
    file = procdir(iJR904, "base_models.bson")
    models = ldat(file; verbose = false);
    model = models[modelkey]
    return uncompress ? ChU.uncompressed_model(model) : model
end

function load_model(modelkey::String, exp::Int; uncompress = true)
    file = procdir(iJR904, "base_models.bson")
    models = ldat(file; verbose = false);
    model = models[modelkey][exp]
    return uncompress ? ChU.uncompressed_model(model) : model
end


text:
function power(maxit::Int64, fidMtx::FidMtx, w::Array{Float64,1}, irz::Array{Int64,1}, irx::Array{Int64,1}, tmax::Float64)
    nz = fidMtx.nz
    nx = fidMtx.nx
    x = randn(nz, nx)
    lambda = 0.0
    for k = 1 : maxit
        b = MultiStepForward(irz, irx, w, x, fidMtx, tmax=tmax)
        y = MultiStepAdjoint(w, b, fidMtx)
        n = norm(vec(x))
        x = y / n
        lambda = n
        println("iteration: $k, maximum eig: $lambda")
    end
    return lambda
end

function fista(shot::Shot, w::Array{Float64,1}, mu::Float64, lambda::Float64, fidMtx::FidMtx, irz::Array{Int64,1}, irx::Array{Int64,1}, tmax::Float64, maxit::Int64)
    nz = fidMtx.nz
    nx = fidMtx.nx
    x = zeros(nz, nx)
    T = mu / (2*lambda)
    t = 1
    yk = copy(x)
    for k = 1 : maxit
        tmpx = copy(x)
        bshot = MultiStepForward(irz, irx, w, yk, fidMtx, tmax=tmax)
        bshot.p = bshot.p - shot.p
        x = MultiStepAdjoint(w, bshot, fidMtx)
        x = yk - x/lambda
        x = vec(x)
        x = softThresh(x, T)
        x = reshape(x, nz, nx)
        tmpt = t
        t = (1 + sqrt(1+4*t^2)) / 2
        yk = x + (tmpt-1)/t * (x-tmpx)
        println("iteration $k")
    end
    return x
end

function softThresh(x::Array{Float64,1}, t::Float64)
    tmp = abs(x) - t
    tmp = (tmp + abs(tmp)) / 2
    y   = sign(x) .* tmp
    return y
end


text:
#Axes
#xlab!/ylab!
for fun in [:x, :y]

    fname = Symbol("$(fun)lab!")

        @eval begin function ($fname)(v::VegaVisualization; title::AbstractString = "", grid::Bool = false, ticks::Number = 0,
                               format::AbstractString = "", formatType::AbstractString="", layer::AbstractString = "front",
                               properties::Dict = Dict{Any, Any}(), tickSize::Number = 0, tickSizeMajor::Number = 0,
                               tickSizeMinor::Number = 0, tickSizeEnd::Number = 0)

            #Find where axis is in array
            a = v.axes[findfirst([z.name == $"$fun" for z in v.scales])]

            #Assign properties, checking for values already set
            title == "" && a.title != "" ? a.title = a.title: a.title = title
            a.grid = grid
            ticks == 0 && a.ticks != nothing? a.ticks = a.ticks: a.ticks = ticks
            tickSize == 0 && a.tickSize != nothing? a.tickSize = a.tickSize: a.tickSize = tickSize
            tickSizeMajor == 0 && a.tickSizeMajor != nothing? a.tickSizeMajor = a.tickSizeMajor: a.tickSizeMajor = tickSizeMajor
            tickSizeMinor == 0 && a.tickSizeMinor != nothing? a.tickSizeMinor = a.tickSizeMinor: a.tickSizeMinor = tickSizeMinor
            tickSizeEnd == 0 && a.tickSizeEnd != nothing? a.tickSizeEnd = a.tickSizeEnd: a.tickSizeEnd = tickSizeEnd
            a.format = format
            a.formatType = formatType
            a.layer = layer

            #See if properties are already defined, and if so merge. If not, define it from function argument
            a.properties == nothing ? properties : merge!(a.properties, properties)

            return v
        end
    end
end

max_stars_repo_name: jpthor/vega


text:
# ok so get some data and import these packages

using BenchmarkTools
using Profile

sample_data=CSV.read("C://Users//tbunc//OneDrive//Data//sample_sequence_data.csv"; copycols=true)

sample_data[:sequence_order] = [1:size(sample_data,1);]
gd = groupby(sample_data, :trace_id)
for subdf in gd
    subdf[:sequence_order,] = [1:size(subdf,1);]
       end

sample_data2 = DataFrame(gd)

Profile.clear()

@btime lev = basic_lev(sample_data2, sequence_id=:trace_id, sequence_event=:activities);
# Before explicit Array declaration:    10.100 ms (79181 allocations: 10.45 MiB)
# After explicit Array declaration:     2.282 ms (12253 allocations: 1.97 MiB)

@btime lev = basic_lev_multi(sample_data2, sequence_id=:trace_id, sequence_event=:activities);
# 15.965 ms (147784 allocations: 14.02 MiB)
x[1]

Profile.print()

# open("C://Users//tbunc//OneDrive//Data//profile.txt", "w") do s
#     Profile.print(IOContext(s, :displaysize => (24, 5000)))
# end


text:
#plotmaker.jl
#I'm doing this in julia because it's easier than in Python
using PGFPlots, DataFrames, CSV, RollingFunctions, LinearAlgebra

df = DataFrame(CSV.File("current_data/dqn.csv"))
name_array = names(df)
# temp = Matrix(df)
# samples = I*temp'
rollouts = Matrix(df)
episodes = (rollouts[:,1]*I)[:]
num_episodes = size(episodes)[1]
reward = (rollouts[:,4]*I)[:] #col 4 for dqn
# reward_std = (rollouts[:,5]*I)[:]
# threshold =
moving_avg = 50
smoothed_reward = rolling(sum, reward, moving_avg)./moving_avg

p = Axis([
    Plots.Scatter(episodes, reward, legendentry="Rollout Data", markSize=0.1),
    Plots.Linear(episodes[moving_avg:end], smoothed_reward, legendentry="Average Reward", mark="none"),
    Plots.Linear(episodes, 300*ones(num_episodes), legendentry="Reward Threshold", mark="none", style="black,dashed")
    ], xlabel="Episode", ylabel="Reward", title="Accumulated Reward vs Episode", xmin=0, xmax=num_episodes
    )
# savefig("plots/biped_ars001.pdf")
# a = Axis(Plots.Linear(x, y, legendentry="My Plot"), xlabel="X", ylabel="Y", title="My Title")
# p = Axis(Plots.Linear(episodes, rollouts, legendentry="reward"))#, markSize=0.1, onlyMarks=true)
# p = Axis(Plots.Linear(episodes, rollouts, legendentry="reward", markSize=0.1, onlyMarks=true), xlabel="Episode", ylabel="Reward", title="Episode Reward vs Episode Number")
p.legendStyle = "at={(1.05,1.0)}, anchor=north west"
save("plots/dqn.pdf", p)


text:
function test_multipermutation_permutations()
    @test collect(permutations(Int[], Int[])) == Array{Int,1}[[]]
    @test length(permutations(Int[], Int[])) == 1
    @test eltype(permutations(Int[], Int[])) == Array{Int,1}

    @test collect(permutations(["a", "b"], [2, 1])) == Array[
        ["a", "a", "b"],
        ["a", "b", "a"],
        ["b", "a", "a"],
    ]
    @test length(permutations(["a", "b"], [2, 1])) == 3
    @test eltype(permutations(["a", "b"], [2, 1])) == Array{String,1}

    @test length(permutations(["M","I","S","P"], [1,4,4,2])) == 34650
    @test length(permutations([:right,:down], BigInt[20,20])) == 137846528820

    @test collect(permutations(Int[], Int[], 0)) == Array{Int,1}[[]]
    @test length(permutations(Int[], Int[], 0)) == 1
    @test eltype(permutations(Int[], Int[], 0)) == Array{Int,1}

    @test collect(permutations(Int[], Int[], 1)) == Array{Int,1}[]
    @test length(permutations(Int[], Int[], 1)) == 0
    @test eltype(permutations(Int[], Int[], 1)) == Array{Int,1}

    @test collect(permutations(["a", "b"], [3, 2], 3)) == Array[
        ["a", "a", "a"],
        ["a", "a", "b"],
        ["a", "b", "a"],
        ["a", "b", "b"],
        ["b", "a", "a"],
        ["b", "a", "b"],
        ["b", "b", "a"],
    ]
    @test length(permutations(["a", "b"], [3, 2], 3)) == 7
    @test eltype(permutations(["a", "b"], [3, 2], 3)) == Array{String,1}

    @test length(permutations(["M","I","S","P"], [1,4,4,2], 11)) == 34650
end

function test_multipermutation_all()
    print(rpad("Math.Combinatorics.Permutation[Multiset]...", 50, ' '))

    test_multipermutation_permutations()

    println("PASS")
end


text:
using CUDA
using SparseArrays
using LinearMaps
#using IterativeSolvers
using Parameters
using BenchmarkTools
using LinearAlgebra
using SBP_GPU
using Random

include("deriv_ops_beta.jl")
include("deriv_ops_GPU.jl")

struct variables
    Nx::Int64
    Ny::Int64
    N::Int64
    hx::Float64
    hy::Float64
    x
    y
    alpha1::Float64
    alpha2::Float64
    alpha3::Float64
    alpha4::Float64
    beta::Float64
end

Nx = 1001
Ny = 1001
N = Nx * Ny

variable = variables(Nx,Ny,N,(1/(Nx-1)),(1/(Ny-1)),0:1/(Nx-1):1,0:1/(Ny-1):1,-1.0,-1.0,-13/(1/(Nx-1)),-13/(1/(Ny-1)),1.0)

# intermediate results
mutable struct intermediates
    Nx::Int64
    Ny::Int64
    N::Int64
    Vol2Face::Array{Array{Float64,1},1}
    du_x::Array{Float64,1}
    du_y::Array{Float64,1}
    du_ops::Array{Float64,1}
    du1::Array{Float64,1}
    du2::Array{Float64,1}
    du3::Array{Float64,1}
    du4::Array{Float64,1}
    du5::Array{Float64,1}
    du6::Array{Float64,1}
    du7::Array{Float64,1}
    du8::Array{Float64,1}
    du9::Array{Float64,1}
    du10::Array{Float64,1}
    du11::Array{Float64,1}
    du12::Array{Float64,1}
    du13::Array{Float64,1}
    du14::Array{Float64,1}
    du15::Array{Float64,1}
    du16::Array{Float64,1}
    du17::Array{Float64,1}
    du0::Array{Float64,1}
    du::Array{Float64,1}
end

intermediate = intermediates(Nx,Ny,N,[zeros(N),zeros(N),zeros(N),zeros(N)],zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N),zeros(N))

mutable struct intermediates_GPU_mutable
    Nx::Int64
    Ny::Int64
    N::Int64
    # du_ops::CuArray{Float64,1}
    du_x::CuArray{Float64,1}
    du_y::CuArray{Float64,1}
    du_ops::CuArray{Float64,1}
    du1::CuArray{Float64,1}
    du2::CuArray{Float64,1}
    du3::CuArray{Float64,1}
    du4::CuArray{Float64,1}
    du5::CuArray{Float64,1}
    du6::CuArray{Float64,1}
    du7::CuArray{Float64,1}
    du8::CuArray{Float64,1}
    du9::CuArray{Float64,1}
    du10::CuArray{Float64,1}
    du11::CuArray{Float64,1}
    du12::CuArray{Float64,1}
    du13::CuArray{Float64,1}
    du14::CuArray{Float64,1}
    du15::CuArray{Float64,1}
    du16::CuArray{Float64,1}
    du17::CuArray{Float64,1}
    du0::CuArray{Float64,1}
    du::CuArray{Float64,1}
end

# N = Nx*Ny
cu_zeros = CuArray(zeros(N))
iGm = intermediates_GPU_mutable(Nx,Ny,N,CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)));

# function myMAT_beta_GPU!(du_GPU::AbstractVector, u_GPU::AbstractVector, container, var_test) # , intermediates_GPU_mutable)
function myMAT_beta_GPU!(du_GPU::AbstractVector, u_GPU::AbstractVector,variable , iGm) # , intermediates_GPU_mutable)
    # @unpack N, y_D2x, y_D2y, y_Dx, y_Dy, y_Hxinv, y_Hyinv, yv2f1, yv2f2, yv2f3, yv2f4, yv2fs, yf2v1, yf2v2, yf2v3, yf2v4, yf2vs, y_Bx, y_By, y_BxSx, y_BySy, y_BxSx_tran, y_BySy_tran, y_Hx, y_Hy = container
    # @unpack Nx,Ny,N,hx,hy,alpha1,alpha2,alpha3,alpha4,beta = var

    # N = Nx*Ny
    # cu_zeros = CuArray(zeros(N))
    # iGm = intermediates_GPU_mutable(Nx,Ny,N,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros,cu_zeros);

    TILE_DIM_1 = 4
    TILE_DIM_2 = 16

    griddim_x = (div(Nx,TILE_DIM_1) + 1, div(Ny,TILE_DIM_2) + 1)
    griddim_y = (div(Nx,TILE_DIM_2) + 1, div(Ny,TILE_DIM_1) + 1)

    blockdim_x = (TILE_DIM_1,TILE_DIM_2)
    blockdim_y = (TILE_DIM_2,TILE_DIM_1)

    # @show typeof(u)
    # @show typeof(iGm.du_x)
    # @show blockdim_x
    # @show griddim_x
    # @show size(u)
    # @show size(iGm.du_x)
    @cuda threads=blockdim_x blocks=griddim_x D2x_GPU_shared(u_GPU,iGm.du_x, variable.Nx, variable.Ny, variable.hx, Val(TILE_DIM_1), Val(TILE_DIM_2))
    # @show Array(iGm.du_x)
    # output = Array(iGm.du_x)
    # output_GPU = iGm.du_x
    synchronize()
    @cuda threads=blockdim_y blocks=griddim_y D2y_GPU_shared(u_GPU,iGm.du_y, variable.Nx, variable.Ny, variable.hy, Val(TILE_DIM_2), Val(TILE_DIM_1))
    synchronize()
    iGm.du_ops = iGm.du_x + iGm.du_y
    output2 = Array(du_ops)
    @cuda threads=blockdim_y blocks=griddim_y BySy_GPU_shared(u_GPU,iGm.du1, variable.Nx, variable.Ny, variable.hy, Val(TILE_DIM_2), Val(TILE_DIM_1))
    # @show iGm.du_x

    synchronize()
    iGm.du2 .= CuArray(VOLtoFACE_beta!(Array(iGm.du1),1,variable.Nx,variable.Ny,variable.N,yv2fs))
    @cuda threads=blockdim_y blocks=griddim_y Hyinv_GPU_shared(iGm.du2,iGm.du3,Nx,Ny,hy, Val(TILE_DIM_2), Val(TILE_DIM_1))
    synchronize()
    iGm.du3 = alpha1 * iGm.du3

    iGm.du5 = VOLtoFACE_beta!(Array(iGm.du1),2,Nx,Ny,N,yv2fs)
    @cuda threads=blockdim_y blocks=griddim_y Hyinv_GPU_shared(iGm.du5,iGm.du6,Nx,Ny,hy, Val(TILE_DIM_2), Val(TILE_DIM_1))
    synchronize()
    iGm.du6 = alpha2 * iGm.du6

    iGm.du7 = CuArray(VOLtoFACE_beta!(Array(u_GPU),3,Nx,Ny,N,yv2fs))
    @cuda threads=blockdim_x blocks=griddim_x BxSx_tran_GPU_shared(iGm.du7,iGm.du8,Nx,Ny,hx,Val(TILE_DIM_1), Val(TILE_DIM_2))
    synchronize()
    @cuda threads=blockdim_x blocks=griddim_x Hxinv_GPU_shared(iGm.du8,iGm.du9,Nx,Ny,hx, Val(TILE_DIM_1), Val(TILE_DIM_2))
    synchronize()
    iGm.du9 = beta * iGm.du9

    @cuda threads=blockdim_x blocks=griddim_x Hxinv_GPU_shared(iGm.du7,iGm.du11,Nx,Ny,hx, Val(TILE_DIM_1), Val(TILE_DIM_2))
    synchronize()
    iGm.du11 =alpha3 * iGm.du11

    iGm.du12 = CuArray(VOLtoFACE_beta!(Array(u_GPU),4,Nx,Ny,N,yv2fs))
    synchronize()
    @cuda threads=blockdim_x blocks=griddim_x BxSx_tran_GPU_shared(iGm.du12,iGm.du13,Nx,Ny,hx,Val(TILE_DIM_1), Val(TILE_DIM_2))
    synchronize()
    @cuda threads=blockdim_x blocks=griddim_x Hxinv_GPU_shared(iGm.du13,iGm.du14,Nx,Ny,hx,Val(TILE_DIM_1), Val(TILE_DIM_2))
    synchronize()
    iGm.du14 = beta * iGm.du14
    @cuda threads=blockdim_x blocks=griddim_x Hxinv_GPU_shared(iGm.du12,iGm.du16,Nx,Ny,hx,Val(TILE_DIM_1), Val(TILE_DIM_2))
    synchronize()
    iGm.du16 = alpha4 * iGm.du16
    synchronize()
    iGm.du0 = iGm.du_ops + iGm.du3 + iGm.du6 + iGm.du9 + iGm.du11 + iGm.du14 + iGm.du16
    synchronize()
    # comment: starting this line, iGm.du17 is not returned with correct solution
    # @cuda threads=blockdim_y blocks=griddim_y Hy_GPU_shared(iGm.du0,iGm.du17,Nx,Ny,hx,Val(TILE_DIM_1),Val(TILE_DIM_2))
    # synchronize()
    @cuda threads=blockdim_x blocks=griddim_x Hy_GPU_shared(iGm.du0,iGm.du17,Nx,Ny,hx,Val(TILE_DIM_1),Val(TILE_DIM_2)) # Looks likes something wrong with griddim_y and blockdim_y
    synchronize()
    @cuda threads=blockdim_x blocks=griddim_x Hx_GPU_shared(iGm.du17,iGm.du,Nx,Ny,hx,Val(TILE_DIM_1),Val(TILE_DIM_2))
    synchronize()
    iGm.du = -1.0 * iGm.du
    # return Array(iGm.du_x)
    # @show output
    # output_final = copy(iGm.du);
    # @show output_final[1:10]
    return Array(iGm.du)
    # return output_final
    # return output2
end

Random.seed!(1234)
u = randn(N);
u_GPU = CuArray(u);

du = similar(u);
du_GPU = similar(u_GPU);

# myMAT_beta_GPU!(du_GPU,u_GPU,iGm,var)

function myMAT_beta!(du::AbstractVector, u::AbstractVector,variable,intermediate)
    D2x_beta!(u,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du_x)
    D2y_beta!(u,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du_y)

    intermediate.du_ops .= intermediate.du_x .+ intermediate.du_y
    BySy_beta!(u,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du1)
    VOLtoFACE_beta!(intermediate.du1,1,variable.Nx,variable.Ny,variable.N,intermediate.Vol2Face)
    Hyinv_beta!(intermediate.Vol2Face[1],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,variable.alpha1,intermediate.du3)     #compute action of P1  .= for faster assignment

    VOLtoFACE_beta!(intermediate.du1,2,variable.Nx,variable.Ny,variable.N,intermediate.Vol2Face)
    Hyinv_beta!(intermediate.Vol2Face[2],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy, variable.alpha2,intermediate.du6)

    #compute action of P2
    VOLtoFACE_beta!(u,3,variable.Nx,variable.Ny,variable.N,intermediate.Vol2Face)
    BxSx_tran_beta!(intermediate.Vol2Face[3],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du8)
    Hxinv_beta!(intermediate.du8,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy, variable.beta, intermediate.du9)
    Hxinv_beta!(intermediate.Vol2Face[3],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,variable.alpha3,intermediate.du11)   #compute action of P3

    VOLtoFACE_beta!(u,4,variable.Nx,variable.Ny,variable.N,intermediate.Vol2Face)
    BxSx_tran_beta!(intermediate.Vol2Face[4],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du13)
    Hxinv_beta!(intermediate.du13,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,variable.beta,intermediate.du14)
    Hxinv_beta!(intermediate.Vol2Face[4],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,variable.alpha4,intermediate.du16)  #compute action of P4
    intermediate.du0 .= intermediate.du_ops .+ intermediate.du3 .+ intermediate.du6 .+ intermediate.du9 .+ intermediate.du11 .+ intermediate.du14 .+ intermediate.du16 #Collect together
    Hy_beta!(intermediate.du0,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du17)
	Hx_beta!(intermediate.du17,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy, -1.0,intermediate.du)
    return intermediate.du
end

function myMAT_beta_no_return!(du::AbstractVector, u::AbstractVector,variable,intermediate)
    D2x_beta!(u,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du_x)
    D2y_beta!(u,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du_y)

    intermediate.du_ops .= intermediate.du_x .+ intermediate.du_y
    BySy_beta!(u,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du1)
    VOLtoFACE_beta!(intermediate.du1,1,variable.Nx,variable.Ny,variable.N,intermediate.Vol2Face)
    Hyinv_beta!(intermediate.Vol2Face[1],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,variable.alpha1,intermediate.du3)     #compute action of P1  .= for faster assignment

    VOLtoFACE_beta!(intermediate.du1,2,variable.Nx,variable.Ny,variable.N,intermediate.Vol2Face)
    Hyinv_beta!(intermediate.Vol2Face[2],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy, variable.alpha2,intermediate.du6)

    #compute action of P2
    VOLtoFACE_beta!(u,3,variable.Nx,variable.Ny,variable.N,intermediate.Vol2Face)
    BxSx_tran_beta!(intermediate.Vol2Face[3],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du8)
    Hxinv_beta!(intermediate.du8,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy, variable.beta, intermediate.du9)
    Hxinv_beta!(intermediate.Vol2Face[3],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,variable.alpha3,intermediate.du11)   #compute action of P3

    VOLtoFACE_beta!(u,4,variable.Nx,variable.Ny,variable.N,intermediate.Vol2Face)
    BxSx_tran_beta!(intermediate.Vol2Face[4],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du13)
    Hxinv_beta!(intermediate.du13,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,variable.beta,intermediate.du14)
    Hxinv_beta!(intermediate.Vol2Face[4],variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,variable.alpha4,intermediate.du16)  #compute action of P4
    intermediate.du0 .= intermediate.du_ops .+ intermediate.du3 .+ intermediate.du6 .+ intermediate.du9 .+ intermediate.du11 .+ intermediate.du14 .+ intermediate.du16 #Collect together
    Hy_beta!(intermediate.du0,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy,intermediate.du17)
	Hx_beta!(intermediate.du17,variable.Nx,variable.Ny,variable.N,variable.hx,variable.hy, -1.0,intermediate.du)
end

function Generate(variable)
    # @unpack Nx,Ny,N,hx,hy,x,y,alpha1,alpha2,alpha3,alpha4,beta = var
    Nx = variable.Nx
    Ny = variable.Ny
    N = variable.N
    hx = variable.hx
    hy = variable.hy
    x = variable.x
    y = variable.y
    alpha1 = variable.alpha1
    alpha2 = variable.alpha2
    alpha3 = variable.alpha3
    alpha4 = variable.alpha4
    beta = variable.beta

    g1 = -pi .* cos.(pi .* x)
    g2 = pi .* cos.(pi .* x .+ pi)
    g3 = sin.(pi .* y)
    g4 = sin.(pi .+ pi .* y)

    f = spzeros(Nx,Ny)
    exactU = spzeros(Nx,Ny)

    for i = 1:Nx
    	for j = 1:Ny
    		f[j,i] = -pi^2 .* sin.(pi .* x[i] + pi .* y[j]) - pi^2 .* sin.(pi .* x[i] + pi .* y[j])
    		exactU[j,i] = sin.(pi .* x[i] + pi .* y[j])
    	end
    end

    f = f[:]
    exact = exactU[:]

    #Construct vector b
    b0 = FACEtoVOL(g1,1,Nx,Ny)
    b1 = alpha1*Hyinv(b0,Nx,Ny,hy)

    b2 = FACEtoVOL(g2,2,Nx,Ny)
    b3 = alpha2*Hyinv(b2,Nx,Ny,hy)

    b4 = FACEtoVOL(g3,3,Nx,Ny)
    b5 = alpha3*Hxinv(b4,Nx,Ny,hx)
    b6 = BxSx_tran(b4,Nx,Ny,hx)
    b7 = beta*Hxinv(b6,Nx,Ny,hx)

    b8 = FACEtoVOL(g4,4,Nx,Ny)
    b9 = alpha4*Hxinv(b8,Nx,Ny,hx)
    b10 = BxSx_tran(b8,Nx,Ny,hx)
    b11 = beta*Hxinv(b10,Nx,Ny,hx)

    bb = b1  + b3  + b5 + b7 + b9 + b11 + f

    #Modify b for PD system
    b12 = Hx(bb,Nx,Ny,hx)
    b = -Hy(b12,Nx,Ny,hy)
    return b,exact
end

b,exact = Generate(variable)

@with_kw struct var_cgs
    Nx = 1001
    Ny = 1001
    N = Nx*Ny
    u = randn(N)
    du = similar(u)
end

var_cg = var_cgs()

#@unpack u, du, N = var_cg
Random.seed!(1234)
du = zeros(N);
u = zeros(N);

# function conjugate(myMAT_beta!,b,container,var,intermediate,maxIteration)
#     @unpack N, y_D2x, y_D2y, y_Dx, y_Dy, y_Hxinv, y_Hyinv, yv2f1, yv2f2, yv2f3, yv2f4, yv2fs, yf2v1, yf2v2, yf2v3, yf2v4, yf2vs, y_Bx, y_By, y_BxSx, y_BySy, y_BxSx_tran, y_BySy_tran, y_Hx, y_Hy = container
#     @unpack Nx,Ny,N,hx,hy,alpha1,alpha2,alpha3,alpha4,beta = var
#     @unpack du_ops,du1,du2,du3,du4,du5,du6,du7,du8,du9,du10,du11,du12,du13,du14,du15,du16,du17,du0 = intermediate

#     u = zeros(N);
#     tol = 1e-16
#     r = b - myMAT_beta!(du,u,container,var,intermediate,maxIteration)
#     p = r
#     rsold = r'*r

#     counts = 0
#     for i = 1:N
#         Ap = myMAT_beta!(du,p,container,var,intermediate,maxIteration)   # can't simply translate MATLAB code, p = r create a link from p to r, once p modified, r will be modified
#         alpha = rsold / (p'*Ap)
#         u = u + alpha * p
#         #axpy!(alpha,p,u)
#         r = r - alpha * Ap
#         #axpy!(-alpha,Ap,r)
#         rsnew = r'*r
#         if sqrt(rsnew) < tol
#             break
#         end
#         p = r + (rsnew/rsold) * p
#         rsold = rsnew;
#         counts += 1
#         #return rsold;
#     end
#     return u, counts
# end

r = similar(u)
function conjugate_beta(myMAT_beta!,r,b,variable,intermediate,maxIteration)
    # @unpack N, y_D2x, y_D2y, y_Dx, y_Dy, y_Hxinv, y_Hyinv, yv2f1, yv2f2, yv2f3, yv2f4, yv2fs, yf2v1, yf2v2, yf2v3, yf2v4, yf2vs, y_Bx, y_By, y_BxSx, y_BySy, y_BxSx_tran, y_BySy_tran, y_Hx, y_Hy = container
    # @unpack Nx,Ny,N,hx,hy,alpha1,alpha2,alpha3,alpha4,beta = var
    # @unpack du_ops,du1,du2,du3,du4,du5,du6,du7,du8,du9,du10,du11,du12,du13,du14,du15,du16,du17,du0 = intermediate

    # N = variable.N
    u = zeros(variable.N);
    du = zeros(variable.N);
    tol = 1e-16

    r .= b .- myMAT_beta!(du,u,variable,intermediate)
    p = copy(r)
    Ap = similar(u)
    rsold = r'*r
    counts = 0
    # maxIteration = 1000
    for i = 1:maxIteration
        Ap .= myMAT_beta!(du,p,variable,intermediate)   # can't simply translate MATLAB code, p = r create a link from p to r, once p modified, r will be modified
        alpha = rsold / (p'*Ap)
        #u = u + alpha * p
        axpy!(alpha,p,u) # BLAS function
        #r = r - alpha * Ap
        axpy!(-alpha,Ap,r)
        rsnew = r'*r
        if sqrt(rsnew) < tol
            break
        end
        #p = r + (rsnew/rsold) * p
        #p .= r .+ (rsnew/rsold) .*p
        p .= (rsnew/rsold) .* p .+ r

        rsold = rsnew;
        counts += 1
        #return rsold;
    end
    return u, counts
end

function conjugate_beta_no_return(myMAT_beta_no_return!,r,b,variable,intermediate,maxIteration)
    # @unpack N, y_D2x, y_D2y, y_Dx, y_Dy, y_Hxinv, y_Hyinv, yv2f1, yv2f2, yv2f3, yv2f4, yv2fs, yf2v1, yf2v2, yf2v3, yf2v4, yf2vs, y_Bx, y_By, y_BxSx, y_BySy, y_BxSx_tran, y_BySy_tran, y_Hx, y_Hy = container
    # @unpack Nx,Ny,N,hx,hy,alpha1,alpha2,alpha3,alpha4,beta = var
    # @unpack du_ops,du1,du2,du3,du4,du5,du6,du7,du8,du9,du10,du11,du12,du13,du14,du15,du16,du17,du0 = intermediate

    # N = variable.N
    u = zeros(variable.N);
    du = zeros(variable.N);
    tol = 1e-16

    # r .= b .- myMAT_beta!(du,u,variable,intermediate)
    myMAT_beta!(du,u,variable,intermediate)
    r .= b .- intermediate.du
    p = copy(r)
    Ap = similar(u)
    rs = [0.0 0.0 0.0]
    # rsold = r'*r
    rs[1] = r'*r
    counts = 0
    # maxIteration = 1000
    for i = 1:maxIteration
        # Ap .= myMAT_beta!(du,p,variable,intermediate)   # can't simply translate MATLAB code, p = r create a link from p to r, once p modified, r will be modified
        myMAT_beta!(du,p,variable,intermediate)   # can't simply translate MATLAB code, p = r create a link from p to r, once p modified, r will be modified
        # alpha = rsold / (p'*intermediate.du)
        alpha = rs[1] / (p'*intermediate.du)
        # rs[3] = rs[1] / (p'*intermediate.du)
        #u = u + alpha * p
        axpy!(alpha,p,u) # BLAS function
        # axpy!(rs[3],p,u) # BLAS function
        #r = r - alpha * Ap
        axpy!(-alpha,intermediate.du,r)
        # axpy!(-rs[3],intermediate.du,r)
        # rsnew = r'*r
        rs[2] = r'*r
        # if sqrt(rsnew) < tol
        if sqrt(rs[2]) < tol
            break
        end
        #p = r + (rsnew/rsold) * p
        #p .= r .+ (rsnew/rsold) .*p
        # p .= (rsnew/rsold) .* p .+ r
        p .= (rs[2]/rs[1]) .* p .+ r
        # axpby!(1., r, rs[2]/rs[1],p)

        # rsold = rsnew;
        rs[1] = rs[2];
        counts += 1
        #return rsold;
    end
    return u, counts
end

# function conjugate_beta_GPU(myMAT_beta_GPU!,r,b,container,var,intermediate,maxIteration)
function conjugate_beta_GPU(myMAT_beta_GPU!,b,var,maxIteration)
    # @unpack N, y_D2x, y_D2y, y_Dx, y_Dy, y_Hxinv, y_Hyinv, yv2f1, yv2f2, yv2f3, yv2f4, yv2fs, yf2v1, yf2v2, yf2v3, yf2v4, yf2vs, y_Bx, y_By, y_BxSx, y_BySy, y_BxSx_tran, y_BySy_tran, y_Hx, y_Hy = container
    # @unpack Nx,Ny,N,hx,hy,alpha1,alpha2,alpha3,alpha4,beta = var
    # @unpack du_ops,du1,du2,du3,du4,du5,du6,du7,du8,du9,du10,du11,du12,du13,du14,du15,du16,du17,du0 = intermediate
    iGm = intermediates_GPU_mutable(Nx,Ny,N,CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)),CuArray(zeros(N)));
    # u = zeros(N);
    # du = zeros(N);
    u_GPU = CuArray(zeros(N))
    du_GPU = CuArray(zeros(N))
    tol = 1e-16

    r = similar(u)
    r .= b .- Array(myMAT_beta_GPU!(du_GPU,u_GPU,iGm,var))
    p = copy(r)
    Ap = similar(u)
    rsold = r'*r
    counts = 0
    # maxIteration = 1000
    for i = 1:maxIteration
        Ap = Array(myMAT_beta_GPU!(du_GPU,CuArray(p),iGm,var))   # can't simply translate MATLAB code, p = r create a link from p to r, once p modified, r will be modified
        Ap = Array(Ap)
        alpha = rsold / (p'*Ap)
        #u = u + alpha * p
        # axpy!(alpha,p,Array(u)) # BLAS function
        u_GPU = u_GPU + alpha * CuArray(p)
        #r = r - alpha * Ap
        # axpy!(-alpha,Ap,r)
        r = r - alpha * Ap
        rsnew = r'*r
        if sqrt(rsnew) < tol
            break
        end
        #p = r + (rsnew/rsold) * p
        #p .= r .+ (rsnew/rsold) .*p
        p .= (rsnew/rsold) .* p .+ r

        rsold = rsnew;
        counts += 1
        #return rsold;
    end
    return u_GPU, counts
end

# (uGPU, countsGPU) = conjugate_beta_GPU(myMAT_beta_GPU!,b,var,100)
# conjugate_beta_GPU(myMAT_beta_GPU!,r,b,container,var,intermediate,100)

(u1,counts1) = conjugate_beta(myMAT_beta!,r,b,variable,intermediate,100)
u1 = copy(u1)
(u1_no_return, counts1_no_return) = conjugate_beta_no_return(myMAT_beta_no_return!,r,b,variable,intermediate,100)
u1_no_return = copy(u1_no_return)
# (u2,counts2) = conjugate_beta(myMAT_beta!,r,b,container,var,intermediate,2000)
# u2 = copy(u2)
# (u3,counts3) = conjugate_beta(myMAT_beta!,r,b,container,var,intermediate,3000)
# u3 = copy(u3)
# (u4,counts4) = conjugate_beta(myMAT_beta!,r,b,container,var,intermediate,1000)
# u4 = copy(u4)
# (u5,counts5) = conjugate_beta(myMAT_beta!,r,b,container,var,intermediate,1000)
# u5 = copy(u5)

counts1
counts2
counts3
counts4
counts5

err_1 = norm(u1 - exact)
err_2 = norm(u2 - exact)
err_3 = norm(u3 - exact)
err_4 = norm(u4 - exact)
err_5 = norm(u5 - exact)

err_norm1 = norm(u1 - exact)
err_norm2 = norm(u2 - exact)

function cg!(du, u, b, myMAT_beta! , tol=1e-16)

    g = similar(u)
    g .= 0
    x = similar(u)
    x .= 0

    myMAT_beta!(g, x, container, var, intermediate)
    g .-= b
    u .=  -g

    gTg = dot(g, g)
    maxiteration = N
    start_time = time()
    for k = 1:maxiteration
        if gTg < tol^2
            end_time = time()
            total_time = end_time - start_time
            return gTg, k, total_time
        end
    gTg = cg_iteration!(x, g, du, u, myMAT_beta!, gTg)
    end

    end_time = time()
    total_time = end_time - start_time
    return u, gTg, maxiteration+1, total_time
end

# function cg_iteration!(u,du, g, w, myMAT_beta!, gTg) version 1
#     myMAT_beta!(du, u, container,var,intermediate)
#     dTw = dot(du, w)
#     alpha = gTg / dTw
#     u .+= alpha .* du
#     g .+= alpha .* w
#     g1Tg1 = dot(g, g)
#     beta = g1Tg1 / gTg
#     u .= .-g .+ beta .* du
#     # u .= u
#     return g1Tg1
# end

function cg_iteration!(x, g, du, u, myMAT_beta!, gTg) #version 2
    #@unpack u, du, N = var_cg
    myMAT_beta!(du, u, container, var, intermediate)
    dTw = dot(u, du)
    alpha = gTg / dTw
    x .+= alpha .* u
    g .+= alpha .* du
    g1Tg1 = dot(g, g)
    beta = g1Tg1 / gTg
    u .= .-g .+ beta .* u
    x .= x
    return g1Tg1
end


text:
using DelayDiffEq

import FiniteDiff
import ForwardDiff

using Test

# Hutchinson's equation
f(u, h, p, t) = p[2] * u * (1 - h(p, t - p[1]) / p[3])
h(p, t) = p[4]

@testset "Gradient" begin
  function test(p)
    prob = DDEProblem(f, p[5], h, eltype(p).((0.0, 10.0)), copy(p);
                      constant_lags = (p[1],))
    sol = solve(prob, MethodOfSteps(Tsit5()); abstol = 1e-14, reltol = 1e-14,
                saveat = 1.0)
    sum(sol)
  end

   # without delay length estimation
  p = [1.5, 1.0, 0.5]
  findiff = FiniteDiff.finite_difference_gradient(p -> test(vcat(1, p, p[end])), p)
  fordiff = ForwardDiff.gradient(p -> test(vcat(1, p, p[end])), p)
  @test maximum(abs.(findiff .- fordiff)) < 1e-6

  # with delay length estimation and without discontinuity
  p = [1.0, 1.5, 1.0, 0.5]
  findiff2 = FiniteDiff.finite_difference_gradient(p -> test(vcat(p, p[end])), p)
  fordiff2 = ForwardDiff.gradient(p -> test(vcat(p, p[end])), p)
  @test_broken maximum(abs.(findiff2 .- fordiff2)) < 2

  # with discontinuity and without delay length estimation
  p = [1.5, 1.0, 0.5, 1.0]
  findiff3 = FiniteDiff.finite_difference_gradient(p -> test(vcat(1, p)), p)
  fordiff3 = ForwardDiff.gradient(p -> test(vcat(1, p)),p)
  @test maximum(abs.(findiff3 .- fordiff3)) < 9

  # consistency checks
  @test findiff2[2:end] ≈ findiff
  @test fordiff2[2:end] ≈ fordiff
  @test_broken findiff3[1:(end-1)] ≈ findiff
  @test_broken fordiff3[1:(end-1)] ≈ fordiff
end

@testset "Jacobian" begin
  function test(p)
    prob = DDEProblem(f, p[5], h, eltype(p).((0.0, 10.0)), copy(p);
                      constant_lags = (p[1],))
    sol = solve(prob, MethodOfSteps(Tsit5()); abstol = 1e-14, reltol = 1e-14,
                saveat = 1.0)
    sol.u
  end

  # without delay length estimation
  p = [1.5, 1.0, 0.5]
  findiff = FiniteDiff.finite_difference_jacobian(p -> test(vcat(1, p, p[end])), p)
  fordiff = ForwardDiff.jacobian(p -> test(vcat(1, p, p[end])), p)
  @test maximum(abs.(findiff .- fordiff)) < 2e-6

  # with delay length estimation and without discontinuity
  p = [1.0, 1.5, 1.0, 0.5]
  findiff2 = FiniteDiff.finite_difference_jacobian(p -> test(vcat(p, p[end])), p)
  fordiff2 = ForwardDiff.jacobian(p -> test(vcat(p, p[end])), p)
  @test_broken maximum(abs.(findiff2 .- fordiff2)) < 3

  # with discontinuity and without delay length estimation
  p = [1.5, 1.0, 0.5, 1.0]
  findiff3 = FiniteDiff.finite_difference_jacobian(p -> test(vcat(1, p)), p)
  fordiff3 = ForwardDiff.jacobian(p -> test(vcat(1, p)),p)
  @test maximum(abs.(findiff3 .- fordiff3)) < 1

  # consistency checks
  @test findiff2[:, 2:end] ≈ findiff
  @test fordiff2[:, 2:end] ≈ fordiff
  @test_broken findiff3[:, 1:(end-1)] ≈ findiff
  @test_broken fordiff3[:, 1:(end-1)] ≈ fordiff
end

@testset "Hessian" begin
  function test(p)
    prob = DDEProblem(f, p[5], h, eltype(p).((0.0, 10.0)), copy(p);
                      constant_lags = (p[1],))
    sol = solve(prob, MethodOfSteps(Tsit5()); abstol = 1e-14, reltol = 1e-14,
                saveat = 1.0)
    sum(sol)
  end

  # without delay length estimation and without discontinuity
  p = [1.5, 1.0, 0.5]
  findiff = FiniteDiff.finite_difference_hessian(p -> test(vcat(1, p, p[end])), p)
  fordiff = ForwardDiff.hessian(p -> test(vcat(1, p, p[end])), p)
  @test maximum(abs.(findiff .- fordiff)) < 1e-5

  # with delay length estimation and without discontinuity
  p = [1.0, 1.5, 1.0, 0.5]
  findiff2 = FiniteDiff.finite_difference_hessian(p -> test(vcat(p, p[end])), p)
  fordiff2 = ForwardDiff.hessian(p -> test(vcat(p, p[end])), p)
  @test_broken maximum(abs.(findiff2 .- fordiff2)) < 25

  # with discontinuity and without delay length estimation
  p = [1.5, 1.0, 0.5, 1.0]
  findiff3 = FiniteDiff.finite_difference_hessian(p -> test(vcat(1, p)), p)
  fordiff3 = ForwardDiff.hessian(p -> test(vcat(1, p)),p)
  @test maximum(abs.(findiff3 .- fordiff3)) < 1

  # consistency checks
  @test findiff2[2:end, 2:end] ≈ findiff
  @test fordiff2[2:end, 2:end] ≈ fordiff
  @test_broken findiff3[1:(end-1), 1:(end-1)] ≈ findiff
  @test_broken fordiff3[1:(end-1), 1:(end-1)] ≈ fordiff
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
using MATLAB
using Test

@test mat"1" == 1
@test mat"[1, 2, 3]" == [1 2 3]

# Test interpolation
x = 1
@test mat"$x + 1" == 2

ret = mat"$y = $x + 2"
@test ret === nothing
@test y == 3

ret = mat"$y = $(x + 3)"
@test ret === nothing
@test y == 4

x = 5
@test mat"$x == 5"

# Test assignment
x = [1, 2, 3, 4, 5]
ret = mat"$x(1:3) = 1"
@test ret === nothing
@test x == [1, 1, 1, 4, 5]
ret = mat"$(x[1:3]) = 2"
@test ret === nothing
@test x == [2, 2, 2, 4, 5]

# Test a more complicated case with assignments on LHS and RHS
x = 20
mat"""
for i = 1:10
   $x = $x + 1;
end
"""

# Test assignment then use
ret = mat"""
$z = 5;
$q = $z;
"""
@test ret === nothing
@test z == 5
@test q == 5

# Test multiple assignment
ret = mat"[$a, $b] = sort([4, 3])"
@test ret === nothing
@test a == [3 4]
@test b == [2 1]

# Test comments
a = 5
@test mat"$a + 1; % = 2" == 6

# Test indexing
c = [1, 2]
@test mat"$c($c == 2)" == 2

# Test line continuations
ret = mat"""
$d ...
= 3
"""
@test ret === nothing
@test d == 3

# Test strings with =
text = "hello = world"
@test mat"strfind($text, 'o = w')" == 5


text:
using PSSFSS
using Test

sh1 = rectstrip(Lx=1, Ly=1.0, Nx=1, Ny=1, Px=1, Py=1, units=inch)

@testset "recstrip" begin
    @test length(sh1.ρ) == 4
    @test sh1.ρ[2] - sh1.ρ[1] == [1,0]
    @test sh1.ρ[3] - sh1.ρ[1] == [0,1]
    @test sh1.ρ[4] - sh1.ρ[1] == [1,1]
    @test sh1.e1 == [1,3,1,2,1]
    @test sh1.e2 == [2,4,3,4,4]
    @test sh1.fe == reshape([2,3,5, 4,5,1], 3,2)
    @test sh1.fv == reshape([1,4,3, 1,2,4], 3,2)
    @test sh1.fufp
    @test sh1.dx == sh1.dy == 0
    @test sh1.rot == 0
    @test sh1.units == inch
    @test sh1.class == 'J'
    @test sh1.s₁ == [1,0]
    @test sh1.s₂ == [0,1]
    @test sh1.β₁ ≈ 2π .* [1,0]
    @test sh1.β₂ ≈ 2π .* [0,1]
end


text:
"""
    Node

Decision tree node type.
"""
mutable struct Node
    value::String
    class::String
    entr::Real
    ig::Real
    children::AbstractArray{Node}
end
Node() = Node("", "", 0, 1, [])
Node(value::String) = Node(value, "", 0, 1, [])
Node(value::String, children::AbstractArray{Node}) =
    Node(value, "", 0, 1, children)

"""Access a node's value."""
value(n::Node) = n.value

"""Access a node's class."""
class(n::Node) = n.class

"""Access a node's entropy."""
entropy(n::Node) = n.entr

"""Access a node's information gain."""
information_gain(n::Node) = n.ig

"""Access a node's children."""
children(n::Node) = n.children

"""Retrieve number of children."""
nchildren(n::Node) = length(n.children)

# Check for equality
Base.:(==)(x::Node, y::Node) =
    (x.value == y.value) & (x.entr == y.entr) & (x.ig == y.ig) &
        (x.children == y.children)

# Add the specified child to the node
Base.push!(n::Node, child::Node) = push!(n.children, child)

# Add all children to the node
Base.append!(n::Node, children::AbstractArray{Node}) =
    append!(n.children, children)

# Display the specified node and its ancestors
Base.show(io::IO, n::Node) = show(io, n)

"""Auxiliary function to display the node and its ancestors with indentation"""
function show(io::IO, n::Node, indent::Integer=0, islast::Bool=false)
    if indent > 1
        print(io, "|   " ^ (indent - 1))
    end
    if indent > 0
        islast ? print(io, "`-- ") : print(io, "|-- ")
    end

    if n.value != ""
        @printf(io, "[%s] (entr=%.4f, ig=%.4f)", n.value, n.entr, n.ig)
    else
        @printf(io, "(%s)", n.class)
    end

    succs = children(n)
    sz = nchildren(n)
    for i in 1:sz
        println(io)
        show(io, succs[i], indent + 1, i == sz)
    end
end


text:
mutable struct robot_state_t <: LCMType
    utime::Int64
    pose::position_3d_t
    twist::twist_t
    num_joints::Int16
    joint_name::Vector{String}
    joint_position::Vector{Float32}
    joint_velocity::Vector{Float32}
    joint_effort::Vector{Float32}
    force_torque::force_torque_t
end

@lcmtypesetup(robot_state_t,
    joint_name => (num_joints,),
    joint_position => (num_joints,),
    joint_velocity => (num_joints,),
    joint_effort => (num_joints,),
)


text:
# tests for deleteHybridizationUpdate! for the 5 taxon network
# we have in the ipad summary
# Claudia September 2014
##########################################################

warn("to run tests_5taxon_delete.jl, you need to set updateGammaz to return true always. this is because we have bad triangles in here")
# types in "types.jl"
include("../src/types.jl")

# needed modules:
using Base.Collections # for updateInCycle with priority queue

# test functions
include("test_functions_5taxon.jl")

tests = ["C","F","G","H","J","D","E","I"];
wrong = String[];
t="C"
for t in tests
    include("../src/functions.jl")
    include("tree_example.jl");
    tp = string("delete_case","$(t).jl");
    println("running $(tp)");
    try
        include(tp)
    catch
        println("error in $(tp)");
        push!(wrong,t);
    end
end

if(!isempty(wrong))
    for t in wrong
        include("../src/functions.jl")
        include("tree_example.jl");
        tp = string("delete_case","$(t).jl");
        println("running $(tp)");
        include(tp)
    end
else
    println("----------NO ERRORS!----------");
end


text:
#!/usr/bin/env julia

using Gtk
using Requests
import Requests: post

const server = "http://localhost"

bld = @GtkBuilder(filename="../global/ui.glade")
win = GAccessor.object(bld, "appWindow")
btn = GAccessor.object(bld, "sendButton")
ibuff = GAccessor.object(bld, "inputBuffer")
obuff = GAccessor.object(bld, "outputBuffer")

signal_connect(win, :destroy) do widget
    Gtk.gtk_quit()
end

function btn_cb(widgetptr::Ptr, data)
    idata = getproperty(ibuff, :text, AbstractString)
    res = post(server,
            headers = Dict("Accept" => "text/plain"),
            data = idata)
    setproperty!(obuff, :text, readall(res))
    nothing
end

signal_connect(btn_cb, btn, "clicked", Void, (), false)

showall(win)
Gtk.gtk_main()

max_stars_repo_name: espectalll/gtk-rest-examples

max_issues_repo_name: espectalll/gtk-rest-examples

max_forks_repo_name: espectalll/gtk-rest-examples


text:
#!/usr/bin/env julia

using WriteVTK

using Random
using Test

const vtk_filename_noext = "polydata"

# Write file with D-dimensional grid.
function write_vtp(Np, D)
    rng = MersenneTwister(42)

    # Workaround changes in rand output in Julia 1.5. See array.jl.
    points = sortslices([rand(rng) for i = 1:D, j = 1:Np], dims=2)

    M = Np >> 2
    @assert 4M == Np

    # Connect every `n` points to form lines.
    lines = let n = 4
        # Try using tuples for connectivity info.
        [MeshCell(PolyData.Lines(), ntuple(d -> i + d, Val(4))) for i = 0:(M - n)]
    end

    # Connect every `n` points to form polygons.
    polys = let n = 5
        [MeshCell(PolyData.Polys(), (i + 1):(i + n)) for i = M:(2M - n)]
    end

    # Each cell is a comination of points (vertices).
    verts = let n = 5
        [MeshCell(PolyData.Verts(), (i + 1):(i + n)) for i = 2M:(3M - n)]
    end

    # Triangle strips
    strips = let n = 12
        [MeshCell(PolyData.Strips(), (i + 1):(i + n)) for i = 3M:(4M - n)]
    end

    # Note: the order of verts, lines, polys and strips is not important.
    # One doesn't even need to pass all of them.
    all_cells = (verts, lines, polys, strips)
    num_cells = sum(length, all_cells)

    fname = vtk_filename_noext * "_$(D)D"
    print("PolyData $(D)D...")
    @time filenames = vtk_grid(fname, points, all_cells...,
                               compress=false, append=false, ascii=true) do vtk
        vtk["my_point_data"] = [randn(rng) for i = 1:3, j = 1:Np]
        vtk["vector_as_tuple"] = ntuple(_ -> [randn(rng) for j = 1:Np], D)
        # NOTE: cell data is not correctly parsed by VTK when multiple kinds of
        # cells are combined in the same dataset.
        # This seems to be a really old VTK issue:
        # - https://vtk.org/pipermail/vtkusers/2004-August/026448.html
        # - https://gitlab.kitware.com/vtk/vtk/-/issues/564
        vtk["my_cell_data"] = [randn(rng) for j = 1:num_cells]
    end
end

function main()
    Np = 128
    filenames = vcat(write_vtp.(Np, 2:3)...)
    println("Saved:  ", join(filenames, "  "))
    filenames
end

main()


text:
input = joinpath(@__DIR__, "input")
lines = readlines(input)

hex = lines[1]
bits = join(map(n -> string(parse(Int, n, base=16), base=2, pad=4), collect(hex)))

function version(bits)
    return parse(Int, bits[1:3], base=2)
end

function header(bits)
    return parse(Int, bits[4:6], base=2)
end

function packets(bits)
    m = parse(Int, bits[7], base=2)
    if m == 1
        return (1, parse(Int, bits[8:18], base=2))
    else
        return (0, parse(Int, bits[8:22], base=2))
    end
end

function val(bits)
    num = bits[7:end]
    acc = []
    i = 1
    while num[i] == '1'
        push!(acc, num[(i + 1):(i + 4)])
        i += 5
    end
    push!(acc, num[(i + 1):(i + 4)])
    rem = num[(i + 5):end]
    return (parse(Int, join(acc), base=2), rem)
end

function mode(bits)
    return parse(Int, bits[7], base=2)
end

function p(bits)
    v = version(bits)
    h = header(bits)
    if h == 4
        n, rem = val(bits)
        return (v, h, n, rem)
    else
        m = mode(bits)
        if m == 1
            acc = []
            _, l = packets(bits)
            rem = bits[19:end]
            while l > 0
                next = p(rem)
                push!(acc, next)
                rem = last(next)
                l -= 1
            end
            return (v, h, acc, rem)
        else
            acc = []
            _, l = packets(bits)
            bs = bits[23:(22 + l)]
            while !isempty(bs)
                next = p(bs)
                push!(acc, next)
                bs = last(next)
            end
            return (v, h, acc, bits[(23 + l):end])
        end
    end
end

pt = p(bits)

function versions(parsetree)
    (v, h, ps, _) = parsetree
    if isempty(ps) || h == 4
        return v
    else
        return v + sum(map(p -> versions(p), ps))
    end
end

function value(parsetree)
    (v, h, ps, _) = parsetree
    if h == 4
        return ps
    elseif h == 0
        return sum(map(value, ps))
    elseif h == 1
        return prod(map(value, ps))
    elseif h == 2
        return minimum(map(value, ps))
    elseif h == 3
        return maximum(map(value, ps))
    elseif h == 5
        vs = map(value, ps)
        vs[1] > vs[2] && return 1
        return 0
    elseif h == 6
        vs = map(value, ps)
        vs[1] < vs[2] && return 1
        return 0
    elseif h == 7
        vs = map(value, ps)
        vs[1] == vs[2] && return 1
        return 0
    end
end

p1 = versions(pt)
p2 = value(pt)

println("-----------------------------------------------------------------------")
println("packet decoder -- part one :: $p1")
println("packet decoder -- part two :: $p2")
println("-----------------------------------------------------------------------")

@assert(p1 == 843)
@assert(p2 == 5390807940351)


text:
using Test, FFTW, Acquisition, GNSSSignals, Random, FFTW
import Unitful: Hz, s

const SAMPLE_CODE = [1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0]

const SAMPLE_CODE2 = [1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, 1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0]

@testset "Estimate signal / noise power" begin
    # Peak in the middle
    power_bins = abs2.(1 / sqrt(2) * complex.(ones(1023,29), ones(1023,29)))
    power_bins[489,11] = power_bins[489,11] + 10^(15 / 10) # 15 dB SNR
    signal_power, noise_power, code_index, doppler_index = Acquisition.est_signal_noise_power(power_bins, 4e6, 1e6)

    @test noise_power ≈ 1
    @test signal_power ≈ 10^(15 / 10)
    @test code_index == 489
    @test doppler_index == 11

    # Peak at doppler and code border
    power_bins = abs2.(1 / sqrt(2) * complex.(ones(1023,29), ones(1023,29)))
    power_bins[1,1] = power_bins[1,1] + 10^(15 / 10) # 15 dB SNR
    signal_power, noise_power, code_index, doppler_index = Acquisition.est_signal_noise_power(power_bins, 4e6, 1e6)

    @test noise_power ≈ 1
    @test signal_power ≈ 10^(15 / 10)
    @test code_index == 1
    @test doppler_index == 1

    power_bins = abs2.(1 / sqrt(2) * complex.(ones(1023,29), ones(1023,29)))
    power_bins[1023,29] = power_bins[1023,29] + 10^(15 / 10) # 15 dB SNR
    signal_power, noise_power, code_index, doppler_index = Acquisition.est_signal_noise_power(power_bins, 4e6, 1e6)

    @test noise_power ≈ 1
    @test signal_power ≈ 10^(15 / 10)
    @test code_index == 1023
    @test doppler_index == 29

    #plotlyjs()
    #surface(-7000:500:7000, 1:1023, power_bins)
    #gui()
end

@testset "Power over code" begin
    integration_time = 2e-3s
    interm_freq = 2e6Hz; code_freq = 1023e3Hz; sample_freq = 4e6Hz; doppler = 1000Hz; f_0 = 1575420e3Hz; code_delay = -5e-4s; range = 1:convert(Int, sample_freq * integration_time)
    carrier =  cis.((2 * π * (interm_freq + doppler) / sample_freq) .* range)
    code_doppler = doppler / (f_0 / code_freq)
    code_phase = -code_delay * (code_freq + code_doppler)
    code_length = length(SAMPLE_CODE)
    code = SAMPLE_CODE[1 .+ mod.(floor.(Int, (code_freq + code_doppler) / sample_freq * range .+ code_phase), code_length)]
    replica_code = SAMPLE_CODE[1 .+ mod.(floor.(Int, code_freq / sample_freq * range), code_length)]
    signal = carrier .* code
    code_freq_dom = fft(replica_code)
    fft_plan = plan_fft(signal)
    power_bins = Acquisition.power_over_code(GPSL1(), signal, fft_plan, code_freq_dom, doppler, sample_freq, interm_freq)
    signal_power, index = findmax(power_bins)
    @test index / sample_freq ≈ code_delay rtol = 2

    #plotlyjs()
    #plot(power_bins)
    #gui()
end

@testset "Power over doppler and code" begin
    integration_time = 2e-3s
    interm_freq = 2e6Hz; code_freq = 1023e3Hz; sample_freq = 4e6Hz; doppler = 1000Hz; f_0 = 1575420e3Hz; code_delay = -5e-4s; range = 1:convert(Int, sample_freq * integration_time)
    carrier =  cis.((2π * (interm_freq + doppler) / sample_freq) .* range)
    code_doppler = doppler / (f_0 / code_freq)
    code_phase = -code_delay * (code_freq + code_doppler)
    code_length = length(SAMPLE_CODE)
    code = SAMPLE_CODE[1 .+ mod.(floor.(Int, (code_freq + code_doppler) / sample_freq * range .+ code_phase), code_length)]
    signal = carrier .* code
    doppler_step = 2 / 3 / integration_time
    doppler_steps = -7000Hz:doppler_step:7000Hz
    power_bins = Acquisition.power_over_doppler_and_code(GPSL1(), signal, 1, doppler_steps, sample_freq, interm_freq)
    signal_power, index = findmax(power_bins)
    c_idx, d_idx = Tuple(CartesianIndices(power_bins)[index])

    @test doppler_step * (d_idx - 1) - 7000Hz == 1000Hz
    @test c_idx / sample_freq ≈ code_delay rtol = 2

    #surface(doppler_steps, 1:4000, power_bins)
    #gui()
    #using PyPlot
    #surf(power_bins)
end

@testset "Acquisition" begin
    Random.seed!(1234)
    integration_time = 2e-3s
    interm_freq = 0.0Hz; code_freq = 1023e3Hz; sample_freq = 4e6Hz; CN0 = 45; doppler = 1000Hz; f_0 = 1575420e3Hz; code_delay = -5e-4s; range = 1:convert(Int, sample_freq * integration_time)
    codes = [SAMPLE_CODE SAMPLE_CODE2]
    #noise_power = 10 * log10(sample_freq)
    noise_power = 1
    signal_power = CN0 - 10 * log10(sample_freq / 1.0Hz)
    noise = 1 / sqrt(2) * complex.(randn(length(range)), randn(length(range)))
    carrier =  cis.((2 * π * (interm_freq + doppler) / sample_freq) .* range)
    code_doppler = doppler / (f_0 / code_freq)
    code_phase = -code_delay * (code_freq + code_doppler)
    code_length = length(SAMPLE_CODE)
    code = SAMPLE_CODE[1 .+ mod.(floor.(Int, (code_freq + code_doppler) / sample_freq * range .+ code_phase), code_length)]
    signal = (carrier .* code) * 10^(signal_power / 20) + noise * 10^(noise_power / 20)
    acq_res = acquire(GPSL1(), signal, sample_freq, 1, dopplers = -7000.0Hz:1 / 3 / (length(signal) / sample_freq):7000.0Hz)
    @test acq_res.carrier_doppler == 1000Hz

    @test acq_res.code_phase ≈ code_phase atol = 1e-3
    @test acq_res.CN0 ≈ CN0 + 10 * log10(integration_time / 1e-3s) atol = 2

    #doppler_step = 2 / 3 / integration_time
    #doppler_steps = -7000:doppler_step:7000
    #using PyPlot
    #surf(acq_res.power_bins)
end


text:
@testset "Range Interface" begin
    @testset "Range Constructors" begin
        @test @inferred(StaticInt(1):StaticInt(10)) == 1:10
        @test @inferred(StaticInt(1):StaticInt(2):StaticInt(10)) == 1:2:10
        @test @inferred(1:StaticInt(2):StaticInt(10)) == 1:2:10
        @test @inferred(StaticInt(1):StaticInt(2):10) == 1:2:10
        @test @inferred(StaticInt(1):2:StaticInt(10)) == 1:2:10
        @test @inferred(1:2:StaticInt(10)) == 1:2:10
        @test @inferred(1:StaticInt(2):10) == 1:2:10
        @test @inferred(StaticInt(1):2:10) == 1:2:10
        @test @inferred(StaticInt(1):UInt(10)) === StaticInt(1):10
        @test @inferred(UInt(1):StaticInt(1):StaticInt(10)) === 1:StaticInt(10)
        @test @inferred(ArrayInterface.OptionallyStaticUnitRange{Int,Int}(1:10)) == 1:10
        @test @inferred(ArrayInterface.OptionallyStaticUnitRange(1:10)) == 1:10

        @inferred(ArrayInterface.OptionallyStaticUnitRange(1:10))

        @test @inferred(ArrayInterface.OptionallyStaticStepRange(StaticInt(1), static(1), static(1))) == 1:1:1
        @test @inferred(ArrayInterface.OptionallyStaticStepRange(StaticInt(1), 1, UInt(10))) == StaticInt(1):1:10
        @test @inferred(ArrayInterface.OptionallyStaticStepRange(UInt(1), 1, StaticInt(10))) == StaticInt(1):1:10
        @test @inferred(ArrayInterface.OptionallyStaticStepRange(1:10)) == 1:1:10

        @test_throws ArgumentError ArrayInterface.OptionallyStaticUnitRange(1:2:10)
        @test_throws ArgumentError ArrayInterface.OptionallyStaticUnitRange{Int,Int}(1:2:10)
        @test_throws ArgumentError ArrayInterface.OptionallyStaticStepRange(1, 0, 10)

        @test @inferred(StaticInt(1):StaticInt(1):StaticInt(10)) === ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10))
        @test @inferred(StaticInt(1):StaticInt(1):10) === ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), 10)
        @test @inferred(1:StaticInt(1):10) === ArrayInterface.OptionallyStaticUnitRange(1, 10)
        @test length(StaticInt{-1}():StaticInt{-1}():StaticInt{-10}()) == 10

        @test UnitRange(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10))) === UnitRange(1, 10)
        @test UnitRange{Int}(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10))) === UnitRange(1, 10)

        @test AbstractUnitRange{Int}(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10))) isa ArrayInterface.OptionallyStaticUnitRange
        @test AbstractUnitRange{UInt}(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10))) isa Base.OneTo
        @test AbstractUnitRange{UInt}(ArrayInterface.OptionallyStaticUnitRange(StaticInt(2), StaticInt(10))) isa UnitRange

        @test @inferred((StaticInt(1):StaticInt(10))[StaticInt(2):StaticInt(3)]) === StaticInt(2):StaticInt(3)
        @test @inferred((StaticInt(1):StaticInt(10))[StaticInt(2):3]) === StaticInt(2):3
        @test @inferred((StaticInt(1):StaticInt(10))[2:3]) === 2:3
        @test @inferred((1:StaticInt(10))[StaticInt(2):StaticInt(3)]) === 2:3

        @test -(StaticInt{1}():StaticInt{10}()) === StaticInt{-1}():StaticInt{-1}():StaticInt{-10}()

        @test reverse(StaticInt{1}():StaticInt{10}()) === StaticInt{10}():StaticInt{-1}():StaticInt{1}()
        @test reverse(StaticInt{1}():StaticInt{2}():StaticInt{9}()) === StaticInt{9}():StaticInt{-2}():StaticInt{1}()
    end

    @test isnothing(@inferred(ArrayInterface.known_first(typeof(1:4))))
    @test isone(@inferred(ArrayInterface.known_first(Base.OneTo(4))))
    @test isone(@inferred(ArrayInterface.known_first(typeof(Base.OneTo(4)))))
    @test isone(@inferred(ArrayInterface.known_first(typeof(StaticInt(1):2:10))))

    @test isnothing(@inferred(ArrayInterface.known_last(1:4)))
    @test isnothing(@inferred(ArrayInterface.known_last(typeof(1:4))))
    @test isone(@inferred(ArrayInterface.known_last(typeof(StaticInt(-1):StaticInt(2):StaticInt(1)))))

    # CartesianIndices
    CI = CartesianIndices((2, 2))
    @test @inferred(ArrayInterface.known_first(typeof(CI))) == CartesianIndex(1, 1)
    @test @inferred(ArrayInterface.known_last(typeof(CI))) == nothing

    CI = CartesianIndices((static(1):static(2), static(1):static(2)))
    @test @inferred(ArrayInterface.known_first(typeof(CI))) == CartesianIndex(1, 1)
    @test @inferred(ArrayInterface.known_last(typeof(CI))) == CartesianIndex(2, 2)

    @test isnothing(@inferred(ArrayInterface.known_step(typeof(1:0.2:4))))
    @test isone(@inferred(ArrayInterface.known_step(1:4)))
    @test isone(@inferred(ArrayInterface.known_step(typeof(1:4))))
    @test isone(@inferred(ArrayInterface.known_step(typeof(Base.Slice(1:4)))))
    @test isone(@inferred(ArrayInterface.known_step(typeof(view(1:4, 1:2)))))

    @testset "length" begin
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(1, 0))) == 0
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(1, 10))) == 10
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), 10))) == 10
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(StaticInt(0), 10))) == 11
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), StaticInt(10)))) == 10
        @test @inferred(length(ArrayInterface.OptionallyStaticUnitRange(StaticInt(0), StaticInt(10)))) == 11

        @test @inferred(length(StaticInt(1):StaticInt(2):StaticInt(0))) == 0
        @test @inferred(length(StaticInt(0):StaticInt(-2):StaticInt(1))) == 0

        @test @inferred(ArrayInterface.known_length(typeof(ArrayInterface.OptionallyStaticStepRange(StaticInt(1), 2, 10)))) === nothing
        @test @inferred(ArrayInterface.known_length(typeof(ArrayInterface.OptionallyStaticStepRange(StaticInt(1), StaticInt(1), StaticInt(10))))) === 10
        @test @inferred(ArrayInterface.known_length(typeof(ArrayInterface.OptionallyStaticStepRange(StaticInt(2), StaticInt(1), StaticInt(10))))) === 9
        @test @inferred(ArrayInterface.known_length(typeof(ArrayInterface.OptionallyStaticStepRange(StaticInt(2), StaticInt(2), StaticInt(10))))) === 5
        @test @inferred(ArrayInterface.known_length(Int)) === 1

        @test @inferred(length(ArrayInterface.OptionallyStaticStepRange(StaticInt(1), 2, 10))) == 5
        @test @inferred(length(ArrayInterface.OptionallyStaticStepRange(StaticInt(1), StaticInt(1), StaticInt(10)))) == 10
        @test @inferred(length(ArrayInterface.OptionallyStaticStepRange(StaticInt(2), StaticInt(1), StaticInt(10)))) == 9
        @test @inferred(length(ArrayInterface.OptionallyStaticStepRange(StaticInt(2), StaticInt(2), StaticInt(10)))) == 5
    end

    @test @inferred(getindex(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), 10), 1)) == 1
    @test @inferred(getindex(ArrayInterface.OptionallyStaticUnitRange(StaticInt(0), 10), 1)) == 0
    @test_throws BoundsError getindex(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), 10), 0)
    @test_throws BoundsError getindex(ArrayInterface.OptionallyStaticStepRange(StaticInt(1), 2, 10), 0)
    @test_throws BoundsError getindex(ArrayInterface.OptionallyStaticUnitRange(StaticInt(1), 10), 11)
    @test_throws BoundsError getindex(ArrayInterface.OptionallyStaticStepRange(StaticInt(1), 2, 10), 11)

    @test ArrayInterface.static_first(Base.OneTo(one(UInt))) === static(1)
    @test ArrayInterface.static_step(Base.OneTo(one(UInt))) === static(1)

    @test Base.setindex(1:5, [6,2], 1:2) == [6,2,3,4,5]
end


text:
module SEED
using Dates, Markdown, Mmap, Printf, SeisIO, SeisIO.FastIO, SeisIO.Formats
path = Base.source_dir()

const id_positions  = Int8[11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
const id_spacer     = 0x2e
const steim         = reverse(collect(0x00000000:0x00000002:0x0000001e), dims=1)
const responses     = Dict{Int64, Any}()
const units_lookup  = Dict{Int64, String}()
const comments      = Dict{Int64, String}()
const abbrev        = Dict{Int64, String}()

# imports
include("SEED/imports.jl")

# files that should be loaded in order
include("SEED/0_seed_read_utils.jl")
include("SEED/1_mSEEDblk.jl")
include("SEED/1_mSEEDdec.jl")
include("SEED/1_dataless_blk.jl")
include("SEED/2_parserec.jl")

# other
include("SEED/dataless.jl")
include("SEED/readmseed.jl")
include("SEED/seed_resp.jl")
include("SEED/seed_support.jl")

# Utils
for i in ls(path*"/SEED/Utils/")
  if endswith(i, ".jl")
    include(i)
  end
end

# exports
include("SEED/exports.jl")

end


text:
using CSV, Distributions, HypothesisTests

data = CSV.read("../data/machine1.csv", header=false)[:,1]
xBar, n = mean(data), length(data)
s = std(data)
alpha = 0.1
t = quantile(TDist(n-1),1-alpha/2)

println("Calculating formula: ", (xBar - t*s/sqrt(n), xBar + t*s/sqrt(n)))
println("Using confint() function: ", confint(OneSampleTTest(xBar,s,n),alpha))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
export VaxInt16, VaxInt32

primitive type VaxInt16 <: VaxInt 16 end
primitive type VaxInt32 <: VaxInt 32 end

VaxInt16(x::UInt16) = reinterpret(VaxInt16, htol(x))
VaxInt16(x::Signed) = reinterpret(VaxInt16, htol(trunc(Int16,x)))

Base.convert(::Type{Int16},x::VaxInt16) = reinterpret(Int16,ltoh(x))
Base.convert(::Type{T},x::VaxInt16) where T <: Union{Int32,Int64,Int128,BigInt,AbstractFloat} = convert(T,convert(Int16,x))

Base.promote_rule(::Type{T},::Type{VaxInt16}) where T <: Union{Int8,VaxInt16} = Int16
Base.promote_rule(::Type{T},::Type{VaxInt16}) where T <: Union{Int16,Int32,Int64,Int128,BigInt,AbstractFloat} = T

Base.promote_type(::Type{VaxInt16}, ::Type{VaxInt16}) = Int16

VaxInt32(x::UInt32) = reinterpret(VaxInt32, htol(x))
VaxInt32(x::Signed) = reinterpret(VaxInt32, htol(trunc(Int32,x)))

Base.convert(::Type{Int32},x::VaxInt32) = reinterpret(Int32,ltoh(x))
Base.convert(::Type{T},x::VaxInt32) where T <: Union{Int16,Int64,Int128,AbstractFloat} = convert(T,convert(Int32,x))

Base.promote_rule(::Type{T},::Type{VaxInt32}) where T <: Union{Int8,Int16,VaxInt16,VaxInt32} = Int32
Base.promote_rule(::Type{T},::Type{VaxInt32}) where T <: Union{Int32,Int64,Int128,BigInt,AbstractFloat} = T

Base.promote_type(::Type{VaxInt32}, ::Type{VaxInt32}) = Int32