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.
#*--------------------------------------------------------------------
#* Author: Tiago Os򱨯rio
#* $Id$
#*-------------------------------------------------------------------
type PIDIncr
	PIDIncr()=begin
		[
			:(PID_output = bias),
			:(diff(PID_dFilt) = 0/"s^2"),
			:(diff(PID_inputFilt)=0/"s"),
			:(diff(PID_setPointFilt)=0/"s"),
		],
		[
			"","","","",
		],
		new(
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Type of PID Incremental",
				:Valid=>["Ideal","Parallel","Series","Ideal_AWBT","Parallel_AWBT","Series_AWBT","Ideal_AW","Parallel_AW","Series_AW"],
				:Default=>"Ideal"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Controller action",
				:Valid=>["Direct","Reverse"],
				:Default=>"Reverse"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Controller mode",
				:Valid=>["Automatic","Manual"],
				:Default=>"Automatic"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Controller mode",
				:Valid=>["Clipped","Unclipped"],
				:Default=>"Clipped"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Derivative term filter constant",
				:Default=>1
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Proportional term setPoint change filter"
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Previous scaled bias",
				:Default=>0.5
			)),
			time_sec(Dict{Symbol,Any}(
				:Brief=>"Derivative time constant"
			)),
			time_sec(Dict{Symbol,Any}(
				:Brief=>"Integral time constant"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Controller gain",
				:Default=>0.5
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Derivative term SP change filter"
			)),
			time_sec(Dict{Symbol,Any}(
				:Brief=>"Input filter time constant"
			)),
			time_sec(Dict{Symbol,Any}(
				:Brief=>"Input filter time constant"
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>1000
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>0
			)),
			control_signal(Dict{Symbol,Any}(
				:Default=>1
			)),
			control_signal(Dict{Symbol,Any}(
				:Protected=>true,
				:PosX=>0,
				:PosY=>0.5
			)),
			control_signal(Dict{Symbol,Any}(
				:Protected=>true,
				:PosX=>0.54,
				:PosY=>1
			)),
			control_signal(),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Derivative term",
				:Unit=>"1/s",
				:Default=>0,
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Derivative term filtered",
				:Default=>0.5,
				:Unit=>"1/s",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Error definition for proportional term",
				:Unit=>"1/s",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Error definition for derivative term",
				:Unit=>"1/s",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Error definition for integral term",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Filtered input",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Integral term",
				:Default=>0,
				:Unit=>"1/s",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Sum of proportional, integral and derivative terms",
				:Unit=>"1/s",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Variable outp scaled between -1 and 1",
				:Hidden=>true,
				:Default=>0.5
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Variable outp",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Proportional term",
				:Default=>0,
				:Unit=>"1/s",
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Filtered setPoint",
				:Default=>0,
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Previous scaled input signal",
				:Default=>0.5,
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Scaled output signal",
				:Default=>0.5,
				:Hidden=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Scaled setPoint",
				:Default=>0.5,
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Integral term multiplier used in anti-reset windup",
				:Hidden=>true
			)),
			DanaReal(Dict{Symbol,Any}(
				:Hidden=>true
			)),
			[
				:(PID_input*(MaxInput - MinInput) = Input - MinInput),
				:(Output = PID_output*(MaxOutput - MinOutput) + MinOutput),
				:(PID_setPoint*(MaxInput - MinInput) = SetPoint - MinInput),
				:((tau + 1e-3*"s")*diff(PID_inputFilt)= PID_input - PID_inputFilt),
				:(tau*diff(PID_inputFilt)= PID_input - PID_inputFilt),
				:((tauSet + 1e-3*"s")*diff(PID_setPointFilt)= PID_setPoint - PID_setPointFilt),
				:(tauSet*diff(PID_setPointFilt)= PID_setPoint - PID_setPointFilt),
				:(PID_error*"s" = PID_inputFilt*(beta-1.0)),
				:(PID_errorD*"s"= PID_inputFilt*(gamma-1.0)),
				:(PID_errorI= 0),
				:(PID_error = beta*diff(PID_setPointFilt) - diff(PID_inputFilt)),
				:(PID_errorD = gamma*diff(PID_setPointFilt) - diff(PID_inputFilt)),
				:(PID_errorI = PID_setPointFilt-PID_inputFilt),
				:(PID_dpropTerm=PID_error),
				:(alpha*(derivTime + 1e-3*"s")*diff(PID_dFilt) = PID_errorD - PID_dFilt),
				:(alpha*(derivTime)*diff(PID_dFilt) = PID_errorD - PID_dFilt),
				:(PID_dderivTerm = derivTime*diff(PID_dFilt)),
				:(diff(PID_outp)=PID_doutp),
				:(PID_outps=2*PID_outp-1),
				:(PID_output=(sign(PID_outps)+1)/2),
				:(PID_output=PID_outp),
				:(PID_output=PID_outp),
				:(PID_action = -1.0),
				:(PID_action = 1.0),
				:(intTime*PID_dintTerm = PID_errorI),
				:(PID_doutp = PID_action*gain*(PID_dpropTerm + PID_dintTerm + PID_dderivTerm)),
				:(PID_AWFactor=1),
				:(intTime*PID_dintTerm = PID_errorI),
				:(PID_doutp = PID_action*(gain*PID_dpropTerm + PID_dintTerm + PID_dderivTerm)),
				:(PID_AWFactor=1),
				:(intTime*PID_dintTerm = PID_errorI),
				:(PID_doutp = PID_action*(gain*(PID_dpropTerm + PID_dintTerm)*(1/"s" + PID_dderivTerm)*"s")),
				:(PID_AWFactor=1),
				:(PID_action*gain*(intTime*PID_dintTerm-PID_errorI) = PID_output-PID_outp),
				:(PID_doutp = PID_action*gain*(PID_dpropTerm + PID_dintTerm + PID_dderivTerm)),
				:(PID_AWFactor=1),
				:(PID_action*gain*(intTime*PID_dintTerm-PID_errorI) = PID_output-PID_outp),
				:(PID_doutp = PID_action*(gain*PID_dpropTerm + PID_dintTerm + PID_dderivTerm)),
				:(PID_AWFactor=1),
				:(PID_action*gain*(intTime*PID_dintTerm-PID_errorI) = PID_output-PID_outp),
				:(PID_doutp = PID_action*(gain*(PID_dpropTerm + PID_dintTerm)*(1/"s" + PID_dderivTerm)*"s")),
				:(PID_AWFactor=1),
				:(intTime*PID_dintTerm = PID_AWFactor*PID_errorI),
				:(PID_doutp = PID_action*gain*(PID_dpropTerm + PID_dintTerm + PID_dderivTerm)),
				:(PID_AWFactor=-tanh(sign(PID_outps)*PID_outps*100-102)),
				:(PID_AWFactor=1),
				:(intTime*PID_dintTerm = PID_AWFactor*PID_errorI),
				:(PID_doutp = PID_action*(gain*PID_dpropTerm + PID_dintTerm + PID_dderivTerm)),
				:(PID_AWFactor=-tanh(sign(PID_outps)*PID_outps*100-102)),
				:(PID_AWFactor=1),
				:(intTime*PID_dintTerm = PID_AWFactor*PID_errorI),
				:(PID_doutp = PID_action*(gain*(PID_dpropTerm + PID_dintTerm)*(1/"s" + PID_dderivTerm)*"s")),
				:(PID_AWFactor=-tanh(sign(PID_outps)*PID_outps*100-102)),
				:(PID_AWFactor=1),
			],
			[
				"Input ","Output ","Set Point ","Input first order filter","Input first order filter","setPoint first order filter","setPoint first order filter","Error definition for proportional term","Error definition for derivative term","Error definition for integral term","Error definition for proportional term","Error definition for derivative term","Error definition for integral term","Calculate proportional term","Derivative term filter","Derivative term filter","Calculate derivative term","Unscaled output","Scale outp","Calculate clipped output when saturated","Calculate clipped output when not saturated","Calculate unclipped output","","","Calculate integral term","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode","Calculate integral term","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode","Calculate integral term","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode","Calculate integral term with anti-windup and bumpless transfer","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode","Calculate integral term with anti-windup and bumpless transfer","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode","Calculate integral term with anti-windup and bumpless transfer","Sum of proportional, integral and derivative terms","Calculate AWFactor - Not in use in this mode","Calculate integral term with anti-windup","Sum of proportional, integral and derivative terms","Calculate AWFactor","Calculate AWFactor","Calculate integral term with anti-windup","Sum of proportional, integral and derivative terms","Calculate AWFactor","Calculate AWFactor","Calculate integral term with anti-windup","Sum of proportional, integral and derivative terms","Calculate AWFactor","Calculate AWFactor",
			],
			[:PID_Select,:Action,:Mode,:Clip,:alpha,:beta,:bias,:derivTime,:intTime,:gain,:gamma,:tau,:tauSet,:MinInput,:MaxInput,:MinOutput,:MaxOutput,],
			[:Input,:Output,:SetPoint,:PID_dderivTerm,:PID_dFilt,:PID_error,:PID_errorD,:PID_errorI,:PID_inputFilt,:PID_dintTerm,:PID_doutp,:PID_outps,:PID_outp,:PID_dpropTerm,:PID_setPointFilt,:PID_input,:PID_output,:PID_setPoint,:PID_AWFactor,:PID_action,]
		)
	end
	PID_Select::DanaSwitcher
	Action::DanaSwitcher
	Mode::DanaSwitcher
	Clip::DanaSwitcher
	alpha::positive
	beta::positive
	bias::control_signal
	derivTime::time_sec
	intTime::time_sec
	gain::positive
	gamma::positive
	tau::time_sec
	tauSet::time_sec
	MinInput::control_signal
	MaxInput::control_signal
	MinOutput::control_signal
	MaxOutput::control_signal
	Input::control_signal
	Output::control_signal
	SetPoint::control_signal
	PID_dderivTerm::control_signal
	PID_dFilt::control_signal
	PID_error::control_signal
	PID_errorD::control_signal
	PID_errorI::control_signal
	PID_inputFilt::control_signal
	PID_dintTerm::control_signal
	PID_doutp::control_signal
	PID_outps::control_signal
	PID_outp::control_signal
	PID_dpropTerm::control_signal
	PID_setPointFilt::control_signal
	PID_input::control_signal
	PID_output::control_signal
	PID_setPoint::control_signal
	PID_AWFactor::DanaReal
	PID_action::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	initials::Array{Expr,1}
	initialNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export PIDIncr
function setEquationFlow(in::PIDIncr)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	if (tau < 1e-3*"s")
		addEquation(4)
	else
		addEquation(5)
	end
	if (tauSet < 1e-3*"s")
		addEquation(6)
	else
		addEquation(7)
	end
	let switch=Mode
		if switch=="Manual"
			addEquation(8)
			addEquation(9)
			addEquation(10)
		elseif switch=="Automatic"
			addEquation(11)
			addEquation(12)
			addEquation(13)
		end
	end
	addEquation(14)
	if (derivTime < 1e-3*"s")
		addEquation(15)
	else
		addEquation(16)
	end
	addEquation(17)
	addEquation(18)
	addEquation(19)
	let switch=Clip
		if switch=="Clipped"
			if abs(PID_outps)>1
				addEquation(20)
			else
				addEquation(21)
			end
			#PID_output = max([0, PID_outp]);

		elseif switch=="Unclipped"
			addEquation(22)
		end
	end
	let switch=Action
		if switch=="Direct"
			addEquation(23)
		elseif switch=="Reverse"
			addEquation(24)
		end
	end
	let switch=PID_Select
		if switch=="Ideal"
			addEquation(25)
			addEquation(26)
			addEquation(27)
		elseif switch=="Parallel"
			addEquation(28)
			addEquation(29)
			addEquation(30)
		elseif switch=="Series"
			addEquation(31)
			addEquation(32)
			addEquation(33)
		elseif switch=="Ideal_AWBT"
			addEquation(34)
			addEquation(35)
			addEquation(36)
		elseif switch=="Parallel_AWBT"
			addEquation(37)
			addEquation(38)
			addEquation(39)
		elseif switch=="Series_AWBT"
			addEquation(40)
			addEquation(41)
			addEquation(42)
		elseif switch=="Ideal_AW"
			addEquation(43)
			addEquation(44)
			if abs(PID_outps)>1 && (PID_action*sign(PID_outps)*PID_errorI)>0
				addEquation(45)
			else
				addEquation(46)
			end
		elseif switch=="Parallel_AW"
			addEquation(47)
			addEquation(48)
			if abs(PID_outps)>1 && (PID_action*sign(PID_outps)*PID_errorI)>0
				addEquation(49)
			else
				addEquation(50)
			end
		elseif switch=="Series_AW"
			addEquation(51)
			addEquation(52)
			if abs(PID_outps)>1 && (PID_action*sign(PID_outps)*PID_errorI)>0
				addEquation(53)
			else
				addEquation(54)
			end
		end
	end
end
function initial(in::PIDIncr)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
end
function atributes(in::PIDIncr,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/PIDincr"
	fields[:Brief]="Model of incremental PIDs."
	fields[:Info]="== Inputs ==
* scaled processs variable.
* scaled bias.
* scaled setpoint.

== Outputs ==
* a scaled output.
"
	drive!(fields,_)
	return fields
end
PIDIncr(_::Dict{Symbol,Any})=begin
	newModel=PIDIncr()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using StatsBase
using StatsBase: PValue, TestStat
using Test, Random, StatsAPI, LinearAlgebra

v1 = [1.45666, -23.14, 1.56734e-13]
v2 = ["Good", "Great", "Bad"]
v3 = [1, 56, 2]
v4 = [-12.56, 0.1326, 2.68e-16]
v5 = [0.12, 0.3467, 1.345e-16]
ct = CoefTable(Any[v1, v2, v3, v4, v5],
               ["Estimate", "Comments", "df", "t", "p"],
               ["x1", "x2", "x3"], 5, 4)
ct_noname = CoefTable(Any[v1, v2, v3, v4, v5],
                      ["Estimate", "Comments", "df", "t", "p"],
                      [], 5, 4)

@test sprint(show, ct) == """
───────────────────────────────────────────────
         Estimate  Comments  df       t       p
───────────────────────────────────────────────
x1    1.45666         Good    1  -12.56  0.1200
x2  -23.14            Great  56    0.13  0.3467
x3    1.56734e-13     Bad     2    0.00  <1e-15
───────────────────────────────────────────────"""

@test sprint(show, ct_noname) == """
────────────────────────────────────────────────
          Estimate  Comments  df       t       p
────────────────────────────────────────────────
[1]    1.45666         Good    1  -12.56  0.1200
[2]  -23.14            Great  56    0.13  0.3467
[3]    1.56734e-13     Bad     2    0.00  <1e-15
────────────────────────────────────────────────"""

@test sprint(show, MIME"text/markdown"(), ct) == """
|    |      Estimate | Comments | df |      t |      p |
|:---|--------------:|---------:|---:|-------:|:-------|
| x1 |   1.45666     |    Good  |  1 | -12.56 | 0.1200 |
| x2 | -23.14        |    Great | 56 |   0.13 | 0.3467 |
| x3 |   1.56734e-13 |    Bad   |  2 |   0.00 | <1e-15 |"""

@test sprint(show, MIME"text/markdown"(), ct_noname) == """
|     |      Estimate | Comments | df |      t |      p |
|:----|--------------:|---------:|---:|-------:|:-------|
| [1] |   1.45666     |    Good  |  1 | -12.56 | 0.1200 |
| [2] | -23.14        |    Great | 56 |   0.13 | 0.3467 |
| [3] |   1.56734e-13 |    Bad   |  2 |   0.00 | <1e-15 |"""

@test length(ct) === 3
@test eltype(ct) ==
    NamedTuple{(:Name, :Estimate, :Comments, :df, :t, :p),
               Tuple{String,Float64,String,Int,Float64,Float64}}
@test collect(ct) == [
    (Name = "x1", Estimate = 1.45666, Comments = "Good", df = 1, t = -12.56, p = 0.12)
    (Name = "x2", Estimate = -23.14, Comments = "Great", df = 56, t = 0.1326, p = 0.3467)
    (Name = "x3", Estimate = 1.56734e-13, Comments = "Bad", df = 2, t = 2.68e-16, p = 1.345e-16)
]

m = [0.11258244478647295  0.05664544616214151  0.38181274408522614  0.8197779704008801
     0.36831406658084287  0.12078054506961555  0.8151038332483567   0.6699313951612162
     0.3444540231363058   0.17957407667101322  0.2422083248151139   0.4530583319523316]
ct = CoefTable(m, ["Estimate", "Stderror", "df", "p"], [], 4)
@test sprint(show, ct) == """
──────────────────────────────────────────
     Estimate   Stderror        df       p
──────────────────────────────────────────
[1]  0.112582  0.0566454  0.381813  0.8198
[2]  0.368314  0.120781   0.815104  0.6699
[3]  0.344454  0.179574   0.242208  0.4531
──────────────────────────────────────────"""
@test length(ct) === 3
@test eltype(ct) ==
    NamedTuple{(:Estimate, :Stderror, :df, :p),
               Tuple{Float64,Float64,Float64,Float64}}
@test collect(ct) == [
    (Estimate = 0.11258244478647295, Stderror = 0.05664544616214151,
     df = 0.38181274408522614, p = 0.8197779704008801)
    (Estimate = 0.36831406658084287, Stderror = 0.12078054506961555,
     df = 0.8151038332483567, p = 0.6699313951612162)
    (Estimate = 0.3444540231363058, Stderror = 0.17957407667101322,
     df = 0.2422083248151139, p = 0.4530583319523316)
]

@test sprint(show, PValue(1.0)) == "1.0000"
@test sprint(show, PValue(1e-1)) == "0.1000"
if VERSION > v"1.6.0-DEV"
    @test sprint(show, PValue(1e-5)) == "<1e-04"
else
    @test sprint(show, PValue(1e-5)) == "<1e-4"
end
@test sprint(show, PValue(NaN)) == "NaN"
@test_throws ErrorException PValue(-0.1)
@test_throws ErrorException PValue(1.1)

@test sprint(show, TestStat(NaN)) == "NaN"
@test sprint(show, TestStat(1e-1)) == "0.10"
@test sprint(show, TestStat(1e-5)) == "0.00"
@test sprint(show, TestStat(π)) == "3.14"

@testset "Union{PValue, TestStat} is Real" begin
    vals = [0.0, Rational(1,3), NaN]
    for T in [PValue, TestStat],
        f in (==, <, ≤, >, ≥, isless, isequal),
        lhs in vals, rhs in vals
        # make sure that T behaves like a Real,
        # regardless of whether it's on the LHS, RHS or both
        @test f(T(lhs), T(rhs)) == f(lhs, rhs)
        @test f(lhs, T(rhs)) == f(lhs, rhs)
        @test f(T(lhs), rhs) == f(lhs, rhs)
    end

    # the (approximate) equality operators get a bit more attention
    for T in [PValue, TestStat]
        @test T(Rational(1,3)) ≈ T(1/3)
        @test Rational(1,3) ≈ T(1/3) atol=0.01
        @test T(Rational(1,3)) isa Real
        @test T(T(0.05)) === T(0.05)
        @test hash(T(0.05)) == hash(0.05)
        @test hash(T(0.05), UInt(42)) == hash(0.05, UInt(42))
    end
end

@test sprint(showerror, ConvergenceException(10)) == "failure to converge after 10 iterations."

@test sprint(showerror, ConvergenceException(10, 0.2, 0.1)) ==
    "failure to converge after 10 iterations. Last change (0.2) was greater than tolerance (0.1)."

@test sprint(showerror, ConvergenceException(10, 0.2, 0.1, "Try changing maxIter.")) ==
    "failure to converge after 10 iterations. Last change (0.2) was greater than tolerance (0.1). Try changing maxIter."

err = @test_throws ArgumentError ConvergenceException(10,.1,.2)
@test err.value.msg == "Change must be greater than tol."

struct MyStatisticalModel <: StatisticalModel
end

StatsAPI.vcov(::MyStatisticalModel) = [1 2; 3 4]
StatsAPI.loglikelihood(::MyStatisticalModel) = 3
StatsAPI.nullloglikelihood(::MyStatisticalModel) = 4
StatsAPI.deviance(::MyStatisticalModel) = 25
StatsAPI.nulldeviance(::MyStatisticalModel) = 40
StatsAPI.dof(::MyStatisticalModel) = 5
StatsAPI.nobs(::MyStatisticalModel) = 100

@testset "StatisticalModel" begin
    m = MyStatisticalModel()

    @test stderror(m) == [1, 2]
    @test aic(m) == 4
    @test aicc(m) ≈ 4.638297872340425
    @test bic(m) ≈ 17.02585092994046
    @test r2(m, :McFadden) ≈ 0.25
    @test r2(m, :CoxSnell) ≈ -0.020201340026755776
    @test r2(m, :Nagelkerke) ≈ 0.24255074155803877
    @test r2(m, :devianceratio) ≈ 0.375

    @test_throws Union{ErrorException, ArgumentError} r2(m, :err)
    @test_throws MethodError r2(m)
    @test adjr2(m, :McFadden) ≈ 1.5
    @test adjr2(m, :devianceratio) ≈ 0.3486842105263158
    @test_throws Union{ErrorException, ArgumentError} adjr2(m, :err)

    @test r2 === r²
    @test adjr2 === adjr²
end

struct MyRegressionModel <: RegressionModel
end

StatsAPI.modelmatrix(::MyRegressionModel) = [1 2; 3 4]

@testset "TestRegressionModel" begin
    m = MyRegressionModel()

    @test crossmodelmatrix(m) == [10 14; 14 20]
    @test crossmodelmatrix(m) isa Symmetric
end

@testset "StatsAPI model reexports" begin
    for f in (fitted, response, responsename, meanresponse,
              modelmatrix, crossmodelmatrix, leverage, cooksdistance, residuals,
              predict, predict!, dof_residual, coef, coefnames, coeftable, confint,
              deviance, islinear, nulldeviance, loglikelihood, nullloglikelihood,
              loglikelihood, loglikelihood, score, nobs, dof, mss, rss,
              informationmatrix, stderror, vcov, weights, isfitted, fit, fit!,
              aic, aicc, bic, r2, r², adjr2, adjr²)
        @test f isa Function
    end
    # Defined but not reexported
    @test StatsBase.params isa Function
    @test StatsBase.params! isa Function
end

ext: jl


text: abstract type Interpolator{T} end


text:
module YaoLang

using LinearAlgebra
using YaoAPI
using YaoCompiler
using IBMQClient
using YaoTargetQASM
using YaoTargetQobj

export @device, @qasm_str, @code_qasm, @code_yao, @ctrl, @measure, @gate

end # module


text:
using MLUtils
using Documenter

# Copy the README to the home page in docs, to avoid duplication.
readme = readlines(joinpath(@__DIR__, "..", "README.md"))
# indexbase = readlines("docs/src/indexbase.md")
open(joinpath(@__DIR__, "src/index.md"), "w") do f
    # for l in indexbase
    #     println(f, l)
    # end
    for l in readme
        println(f, l)
    end
end

DocMeta.setdocmeta!(MLUtils, :DocTestSetup, :(using MLUtils); recursive=true)

makedocs(;
    modules=[MLUtils],
    doctest=true, clean=true,
    sitename = "MLUtils.jl",
    pages = ["Home" => "index.md",
             "API" => "api.md"],
)

deploydocs(repo="github.com/JuliaML/MLUtils.jl.git")


text:
using FDModBiotCyl
using Dierckx

#=========================
Creating a homogeneous model
=========================#
function makemodel_homogeneous_Elastic()
#Necessary input matrices to the main loop function are:
#Rho: Bulk density
#Rhof: Fluid density
#M,C,H: Poroelastic moduli
#G: Formation shear moduli
#D1,D2: Poroelastic moduli in Ou's formulation
#Flag_AC,Flag_E: flag specifying acoustic region or elastic region
# This function creates the above input parameters assuming a
# homogeneous elastic media

#Model size
   nr=251 #samples
   nz=351 #samples
   dr=0.5 #meter
   dz=0.5 #meter

   #==============================================
   Initializing poroelastic parameters (Sidler's)
   ==============================================#
   #First creating Sidler's poroelastic parameters and then converting to
   #necessary parameters for our FD
   #----solid phase-----
   Km=zeros(nz,nr) #frame
   Ks=zeros(nz,nr) #grain
   G=zeros(nz,nr) #bulk
   Rho=zeros(nz,nr) #bulk (weighted average of Rhos and Rhof)
   Rhos=zeros(nz,nr) #grain
   #---fluid phase----
   Kf=zeros(nz,nr)
   Rhof=zeros(nz,nr)
   Kappa0=zeros(nz,nr) #static permeability (m^2, 1D=9.869E-13m^2)
   Eta=zeros(nz,nr) #dynamic viscosity Pa.s=kg/(m.s)
   #---other parameters---
   Phi=zeros(nz,nr)
   Tot=zeros(nz,nr) #Tortuosity factor (see Sidler 2014)
   #End initializing poroelastic parameters

   #==============================================
   Fluid parameters
   ==============================================#
   #Assuming that fluid properties are identical everywhere
   Kf[:,:]=ones(nz,nr)*2.25*10^9
   Rhof[:,:]=ones(nz,nr)*1000
   Eta[:,:]=ones(nz,nr)*1.0*10^(-3) #water 1E-3 Pa.s

   Flag_AC=zeros(nz,nr)

   #==========================
   Homogeneous Elastic model
   ==========================#
   Flag_E=zeros(nz,nr)
   Vp1_elastic=2500.0 #Or, Specify K_elastic
   Vs1_elastic=2000.0
   Rho1_elastic=2000.0

   G_elastic=Vs1_elastic^2*Rho1_elastic
   K_elastic=Vp1_elastic^2*Rho1_elastic-4/3*G_elastic #Activate here if you have specified Vp_elastic

   Flag_E=ones(nz,nr)
   Rho=Rho1_elastic*ones(nz,nr) #Rho=Rho_elastic
   G=ones(nz,nr)*G_elastic #G=Gs (Ou's)
   #--
   Km=ones(nz,nr)*K_elastic #Km
   Ks=ones(nz,nr)*K_elastic #Ks=Km (-> M=Inf, Ou's)
   #--

   #=====================================================
   Converting Sidler's poroelastic parameters into Ou's
   =====================================================#
   #Rho, Rhof, M, C, H, G, D1, D2
   #Rho as is
   #Rhof as is
   Alpha=ones(nz,nr)-Km./Ks
   M=( (Alpha-Phi)./Ks+Phi./Kf ).^(-1)
   C=M.*Alpha
   H=Km+4/3*G+M.*Alpha.^2
   #G as is
   D1=Tot.*Rhof./Phi #(1+2/m)*alpha_inf*rhof/phi
   D2=Eta./Kappa0

   #Correction of H and C for elastic media
   for ir=1:nr
      for iz=1:nz
         if(Flag_E[iz,ir]==1)
            H[iz,ir]=Ks[iz,ir]+4/3*G[iz,ir]
            C[iz,ir]=Ks[iz,ir] #Ou's
         end

         if(Flag_AC[iz,ir]==1)
#            D1[iz,ir]=1.0 #Guan's
            D1[iz,ir]=1000. #Ou's D1=rhof
            D2[iz,ir]=0. #Ou's D2=0
         end
      end
   end

   return nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,Flag_AC,Flag_E
end

#=========================
Green's function
=========================#
function Green_Aki(vp,vs,rho,tvec,src_func,srcgeom,recgeom)
#Aki and Rechards eq 4.23
#vp=2500
#vs=2000
#rho=2500
println("Aki and Rechards Green's function ...")

tmp_R=recgeom-repeat(srcgeom,nrec,1) #(z,r)
R=map(x->sqrt(x),sum(tmp_R.^2,dims=2))
gammax=tmp_R[:,2]./R
gammaz=tmp_R[:,1]./R
src_spl=Spline1D(tvec[:],src_func[:],k=1)
    #================
          Gzz
    ================#
     amp_p=1/(4pi*rho*vp^2)*gammaz.*gammaz./R
     amp_s=-1/(4pi*rho*vs^2)*(gammaz.*gammaz-ones(size(gammax)))./R
     Gzz_far=zeros(nt,nrec)
     for irec=1:nrec
        tmp=amp_p[irec]*src_spl(tvec[:]-R[irec]/vp*ones(size(tvec)))+
            amp_s[irec]*src_spl(tvec[:]-R[irec]/vs*ones(size(tvec)))
        Gzz_far[:,irec]=tmp #far-field terms only (displacement-body force)
     end
     #inclusion of a near-field term
     Gzz_near=zeros(nt,nrec)
     amp_nf=1/(4pi*rho)*(3*gammaz.*gammaz-ones(size(gammax)))./(R.^3)

     for irec=1:nrec
        tau=range(R[irec]/vp,R[irec]/vs,length=nt) # range object (no memory allocation)
        tau=collect(tau) # a vector
        dt_tau=tau[2]-tau[1]
        tmp_nf=zeros(nt,1)
        for it=1:nt
           src_intp=src_spl(tvec[it]*ones(size(tau))-tau[:])
           tmp_nf[it]=sum(src_intp.*tau)*dt_tau
        end
        Gzz_near[:,irec]=tmp_nf*amp_nf[irec]
     end

     Gzz=Gzz_far+Gzz_near

#converting displacement to particle velocity
  tmp=(Gzz[2:end,:]-Gzz[1:end-1,:])/dt
  vz_an=zeros(nt,nrec)
  for irec=1:nrec
    tmp_spl=Spline1D(tvec[1:end-1]+dt/2*ones(size(tmp[:,irec])),tmp[:,irec],k=1)
    vz_an[:,irec]=tmp_spl(tvec)
  end

  return vz_an
end

#----MAIN loop function with src----
function main_loop!(nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,dt,nt,T,
   vr,vz,trr,tpp,tzz,trz,
   vfr,vfz,pf,
   src_index,src_dn,src_func,
   Flag_AC,Flag_E,
   nrec,rec_vr,rec_vz,index_allrec_vr,index_allrec_vz,
   LPML_r,LPML_z,PML_Wr,PML_Wz,PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2,
   snapshots_vr,snapshots_vz,snapshots_trr,nsnap,itvec_snap,nskip)

   println("FD main loop ...")

#--
Flag_vf_zero=get_Flag_vf_zero(Flag_AC,Flag_E,nr,nz)

#initialize matrices of PML
Prz_T,Pzz_T,Rz_T,Srz_T,PzzPE_T,RzPE_T,
Prz_B,Pzz_B,Rz_B,Srz_B,PzzPE_B,RzPE_B,
Prr_R,Qrp_R,Pzr_R,Qzp_R,Rr_R,Rp_R,Rrz_R,PrrPE_R,RrPE_R,RpPE_R,
Prr_TR,Qrp_TR,Prz_TR,Pzr_TR,Qzp_TR,Pzz_TR,Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,PrrPE_TR,PzzPE_TR,RrPE_TR,RpPE_TR,RzPE_TR,
Prr_BR,Qrp_BR,Prz_BR,Pzr_BR,Qzp_BR,Pzz_BR,Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,PrrPE_BR,PzzPE_BR,RrPE_BR,RpPE_BR,RzPE_BR=
init_PML_Por(LPML_r,LPML_z,nr,nz)

memT_vr,memT_vz,memT_trr,memT_tpp,memT_tzz,memT_trz,memT_vfr,memT_vfz,memT_pf,
memB_vr,memB_vz,memB_trr,memB_tpp,memB_tzz,memB_trz,memB_vfr,memB_vfz,memB_pf,
memR_vr,memR_vz,memR_trr,memR_tpp,memR_tzz,memR_trz,memR_vfr,memR_vfz,memR_pf=
init_memory_variables_Por(LPML_r,LPML_z,nr,nz)
#return

#required for BC
vr_old=zeros(nz,nr)
vz_old=zeros(nz,nr)
vfr_old=zeros(nz,nr)
vfz_old=zeros(nz,nr)

for ii=1:nt

#----Time at (ii-1)*dt---(updating velocities)-----

#--PML: save velocity at previous step
PML_save_vel!(memT_vr,memT_vz,memT_vfr,memT_vfz,
              memB_vr,memB_vz,memB_vfr,memB_vfz,
              memR_vr,memR_vz,memR_vfr,memR_vfz,
              vr,vz,vfr,vfz,nr,nz,LPML_z,LPML_r)

# Main velocity update!
update_velocity_1st_Por!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt,LPML_z,LPML_r)

#PML: update velocity
PML_update_vel!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
              Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt,
              Prz_T,Pzz_T,PzzPE_T,
              Prz_B,Pzz_B,PzzPE_B,
              Prr_R,Qrp_R,Pzr_R,Qzp_R,Rr_R,Rp_R,Rrz_R,
              PrrPE_R,RrPE_R,RpPE_R,
              Prz_TR,Pzz_TR,PzzPE_TR,
              Prr_TR,Qrp_TR,Pzr_TR,Qzp_TR,Rr_TR,Rp_TR,Rrz_TR,
              PrrPE_TR,RrPE_TR,RpPE_TR,
              Prz_BR,Pzz_BR,PzzPE_BR,
              Prr_BR,Qrp_BR,Pzr_BR,Qzp_BR,Rr_BR,Rp_BR,Rrz_BR,
              PrrPE_BR,RrPE_BR,RpPE_BR,
              LPML_z,LPML_r)

# Making sure RightBC for velocity (zero values are vr)
ApplyBCRight_vel!(vr,vz,
      trr,tpp,tzz,trz,
      vfr,vfz,pf,
      Rho,Rhof,D1,D2,Flag_vf_zero,nr,nz,dr,dz,dt)

#---velocity B.Cs-----
#Periodic Left Edge
ApplyBCLeft_vel!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                          Rho,Rhof,D1,D2,Flag_vf_zero,
                          nr,nz,dr,dz,dt,
                          Prz_T,Pzz_T,PzzPE_T,
                          Prz_B,Pzz_B,PzzPE_B,
                          LPML_z)

#PML: update memory variables for stress (Rx and Sxx) using velocity at two time steps
PML_update_memRS!(Rz_T,Srz_T,RzPE_T,
                  Rz_B,Srz_B,RzPE_B,
                  Rr_R,Rp_R,Rrz_R,RrPE_R,RpPE_R,
                  Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,RrPE_TR,RpPE_TR,RzPE_TR,
                  Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,RrPE_BR,RpPE_BR,RzPE_BR,
                  memT_vr,memT_vz,memT_vfr,memT_vfz,
                  memB_vr,memB_vz,memB_vfr,memB_vfz,
                  memR_vr,memR_vz,memR_vfr,memR_vfz,
                  vr,vz,vfr,vfz,nr,nz,dr,dz,dt,LPML_z,LPML_r,PML_Wz,PML_Wr,PML_IWr,PML_Wz2,PML_Wr2,PML_IWr2)

#--src injection (velocity)
srcamp=src_func[ii]
srcamp=srcamp*dt/Rho[1,1]
srcapply!(vz,src_index,src_dn,srcamp)

#----Time at (ii-1)*dt+dt/2----(updating stress)---
#--PML: save stress at previous step
PML_save_stress!(memT_trr,memT_tpp,memT_tzz,memT_trz,memT_pf,
                 memB_trr,memB_tpp,memB_tzz,memB_trz,memB_pf,
                 memR_trr,memR_tpp,memR_tzz,memR_trz,memR_pf,
                 trr,tpp,tzz,trz,pf,nr,nz,LPML_z,LPML_r)

#println("update stress")

# Main stress update!
update_stress_1st_Por!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
     M,C,H,G,nr,nz,dr,dz,dt,LPML_z,LPML_r)
#error()

#PML: stress update
PML_update_stress!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                  M,C,H,G,nr,nz,dr,dz,dt,
                  Rz_T,Srz_T,RzPE_T,
                  Rz_B,Srz_B,RzPE_B,
                  Rr_R,Rp_R,Rrz_R,RrPE_R,RpPE_R,
                  Rr_TR,Rp_TR,Rz_TR,Rrz_TR,Srz_TR,RrPE_TR,RpPE_TR,RzPE_TR,
                  Rr_BR,Rp_BR,Rz_BR,Rrz_BR,Srz_BR,RrPE_BR,RpPE_BR,RzPE_BR,
                  LPML_z,LPML_r)

# Making sure RightBC for stress (zero values are trp and trz)
ApplyBCRight_stress!(vr,vz, #Use with flag_zero=1 (see ApplyBCRight_stress1D01)
     trz,
     G,nr,nz,dr,dz,dt)

#---stress B.Cs
#println("update stress Left")
ApplyBCLeft_stress!(vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf,
                    M,C,H,G,nr,nz,dr,dz,dt,
                    Rz_T,Srz_T,RzPE_T,
                    Rz_B,Srz_B,RzPE_B,
                    LPML_z)

#--PML: update memory variables for velocity (Pxx and Qxx) using stress at two time steps
PML_update_memPQ!(Prz_T,Pzz_T,PzzPE_T,
                  Prz_B,Pzz_B,PzzPE_B,
                  Prr_R,Qrp_R,Pzr_R,Qzp_R,PrrPE_R,
                  Prr_TR,Qrp_TR,Prz_TR,Pzr_TR,Qzp_TR,Pzz_TR,PrrPE_TR,PzzPE_TR,
                  Prr_BR,Qrp_BR,Prz_BR,Pzr_BR,Qzp_BR,Pzz_BR,PrrPE_BR,PzzPE_BR,
                  memT_trr,memT_tpp,memT_tzz,memT_trz,memT_pf,
                  memB_trr,memB_tpp,memB_tzz,memB_trz,memB_pf,
                  memR_trr,memR_tpp,memR_tzz,memR_trz,memR_pf,
                  trr,tpp,tzz,trz,pf,nr,nz,dr,dz,dt,LPML_z,LPML_r,PML_Wz,PML_Wr,PML_IWr,PML_Wz2,PML_Wr2,PML_IWr2)

#-----End of updating field----

#println("receiver")
#Receiver field (extraction)
getRecData_from_index!(vr,rec_vr,index_allrec_vr,nrec,ii)
getRecData_from_index!(vz,rec_vz,index_allrec_vz,nrec,ii)

end #i Time

   println("FD main loop ... finished.")

end #function
#------------------------------

#-----ENTRY POINT HERE----
println("===============================")
println("Homogeneous elastic model test.")
println("===============================")
#======================
General variables
======================#

#time samples
dt=1E-4
nt=751
T=(nt-1)*dt
tvec=range(0.0,T,length=nt) # range object (no memory allocation)
tvec=collect(tvec) # a vector

#PML thickness in samples
LPML_r=20
LPML_z=20

#======================
Creating a model
======================#
nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,Flag_AC,Flag_E=makemodel_homogeneous_Elastic()

#======================
Stability check
======================#
Vmax=maximum((H./Rho).^(0.5))
check_stability01(dt,dr,Vmax,0)
check_stability01(dt,dz,Vmax,0)

#======================
Creating a src wavelet
======================#
f0=100 #src Freq
delay=1/f0*1.0
src_func=myricker2(tvec,f0,delay,1) #when using 1st derivative Gaussian
tmp_maxamp=maximum(map(abs,src_func))
src_func=src_func/tmp_maxamp

#================================
Point src geometry
================================#
srcgeom=zeros(1,2) #(z,r)
srcgeom[1,1]=(nz-1)*dz/2 #z meter
srcgeom[1,2]=0*dr #r meter.
#Gaussian point src
wsize=7 #window size (odd number)
wsigma=1 #std
src_index,src_dn=get_srcindex_pGauss(srcgeom,dr,dz,wsize,wsigma)

#==============================
Initializing field variables
==============================#
vr,vz,trr,tpp,tzz,trz,vfr,vfz,pf=init_fields_Por(nz,nr) #
# Initializing PML field variables
LPML_r,LPML_z,PML_Wr,PML_Wz,
PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2=init_PML_profile(LPML_r,LPML_z,Vmax,dr,dz,nr,f0)

#===================
Receiver geometry
====================#
nrec=nz
recgeom=zeros(nrec,2) #(z,r)
for irec=1:nrec
   recgeom[irec,1]=dz*(irec-1)
   recgeom[irec,2]=25.0
end
rec_vr,rec_vz,index_allrec_vr,index_allrec_vz=init_receiver_geophone(recgeom,nrec,dr,dz,nr,nz,nt)

#==============================
Snapshot settings
==============================#
nskip=30 #change here to adjust time sampling (1~nt)
snapshots_vz,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)
snapshots_vr,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)
snapshots_trr,nsnap,itvec_snap=init_snapshots(nt,nz,nr,nskip)

#==============================
Start main FD Loop
==============================#
main_loop!(nr,nz,dr,dz,Rho,Rhof,M,C,H,G,D1,D2,dt,nt,T,
   vr,vz,trr,tpp,tzz,trz,
   vfr,vfz,pf,
   src_index,src_dn,src_func,
   Flag_AC,Flag_E,
   nrec,rec_vr,rec_vz,index_allrec_vr,index_allrec_vz,
   LPML_r,LPML_z,PML_Wr,PML_Wz,PML_IWr,PML_Wr2,PML_Wz2,PML_IWr2,
   snapshots_vr,snapshots_vz,snapshots_trr,nsnap,itvec_snap,nskip)

#Checking with theoretical Green's function
rec_vz_Aki=Green_Aki(2500,2000,2000,tvec,src_func,srcgeom,recgeom) #vp,vs,rho

tol=0.05
@test sqrt(sum((rec_vz[:,151] - rec_vz_Aki[:,151]) .^ 2)) / sqrt(sum(rec_vz[:,151] .^ 2)) < tol


text:
import GPLinearODEMaker.powers_of_negative_one

"""
    se_se_kernel(hyperparameters, δ, dorder; shift_ind=0)

Created by kernel_coder(). Requires 3 hyperparameters.
Likely created using se_se_kernel_base() as an input.
Use with include("src/kernels/se_se_kernel.jl").

# Arguments
- `hyperparameters::Vector`: The hyperparameter values. For this kernel, they should be `["λ1", "λ2", "sratio"]`
- `δ::Real`: The difference between the inputs (e.g. `t1 - t2`)
- `dorder::Vector{<:Integer}`: How many times to differentiate with respect to the inputs and the `hyperparameters` (e.g. `dorder=[0, 1, 0, 2]` would correspond to differentiating once w.r.t the second input and twice w.r.t `hyperparameters[2]`)
- `shift_ind::Integer=0`: If changed, the index of which hyperparameter is the `δ` shifting one
"""
function se_se_kernel(
    hyperparameters::Vector{<:Real},
    δ::Real,
    dorder::Vector{<:Integer};
    shift_ind::Integer=0)

    include_shift = shift_ind!=0
    @assert length(hyperparameters)==3+Int(include_shift) "hyperparameters is the wrong length"
    dorder_len = 5 + Int(include_shift)
    @assert length(dorder)==dorder_len "dorder is the wrong length"
    @assert maximum(dorder) < 3 "No more than two derivatives for each time or hyperparameter can be calculated"

    dorder2 = dorder[2]
    dorder[2] += dorder[1]

    if include_shift
        dorder[2] += dorder[shift_ind+2]
        δ += hyperparameters[shift_ind]
        hyperparameters_view = view(hyperparameters, 1:4 .!= shift_ind)
        dorder_view = view(dorder, 2:dorder_len)
        dorder_view = view(dorder_view, 1:(dorder_len-1) .!= (shift_ind+1))
    else
        hyperparameters_view = hyperparameters
        dorder_view = view(dorder, 2:dorder_len)
    end

    λ1 = hyperparameters_view[1]
    λ2 = hyperparameters_view[2]
    sratio = hyperparameters_view[3]

    if dorder_view==[6, 0, 0, 2]
        func = -30*exp((-1/2)*δ^2/λ2^2)/λ2^6 + 90*exp((-1/2)*δ^2/λ2^2)*δ^2/λ2^8 - 30*exp((-1/2)*δ^2/λ2^2)*δ^4/λ2^10 + 2*exp((-1/2)*δ^2/λ2^2)*δ^6/λ2^12
    end

    if dorder_view==[5, 0, 0, 2]
        func = -30*exp((-1/2)*δ^2/λ2^2)*δ/λ2^6 + 20*exp((-1/2)*δ^2/λ2^2)*δ^3/λ2^8 - 2*exp((-1/2)*δ^2/λ2^2)*δ^5/λ2^10
    end

    if dorder_view==[4, 0, 0, 2]
        func = 6*exp((-1/2)*δ^2/λ2^2)/λ2^4 - 12*exp((-1/2)*δ^2/λ2^2)*δ^2/λ2^6 + 2*exp((-1/2)*δ^2/λ2^2)*δ^4/λ2^8
    end

    if dorder_view==[3, 0, 0, 2]
        func = 6*exp((-1/2)*δ^2/λ2^2)*δ/λ2^4 - 2*exp((-1/2)*δ^2/λ2^2)*δ^3/λ2^6
    end

    if dorder_view==[2, 0, 0, 2]
        func = -2*exp((-1/2)*δ^2/λ2^2)/λ2^2 + 2*exp((-1/2)*δ^2/λ2^2)*δ^2/λ2^4
    end

    if dorder_view==[1, 0, 0, 2]
        func = -2*exp((-1/2)*δ^2/λ2^2)*δ/λ2^2
    end

    if dorder_view==[0, 0, 0, 2]
        func = 2*exp((-1/2)*δ^2/λ2^2)
    end

    if dorder_view==[6, 0, 1, 1]
        func = 180*exp((-1/2)*δ^2/λ2^2)*sratio/λ2^7 - 750*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio/λ2^9 + 390*exp((-1/2)*δ^2/λ2^2)*δ^4*sratio/λ2^11 - 54*exp((-1/2)*δ^2/λ2^2)*δ^6*sratio/λ2^13 + 2*exp((-1/2)*δ^2/λ2^2)*δ^8*sratio/λ2^15
    end

    if dorder_view==[5, 0, 1, 1]
        func = 180*exp((-1/2)*δ^2/λ2^2)*δ*sratio/λ2^7 - 190*exp((-1/2)*δ^2/λ2^2)*δ^3*sratio/λ2^9 + 40*exp((-1/2)*δ^2/λ2^2)*δ^5*sratio/λ2^11 - 2*exp((-1/2)*δ^2/λ2^2)*δ^7*sratio/λ2^13
    end

    if dorder_view==[4, 0, 1, 1]
        func = -24*exp((-1/2)*δ^2/λ2^2)*sratio/λ2^5 + 78*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio/λ2^7 - 28*exp((-1/2)*δ^2/λ2^2)*δ^4*sratio/λ2^9 + 2*exp((-1/2)*δ^2/λ2^2)*δ^6*sratio/λ2^11
    end

    if dorder_view==[3, 0, 1, 1]
        func = -24*exp((-1/2)*δ^2/λ2^2)*δ*sratio/λ2^5 + 18*exp((-1/2)*δ^2/λ2^2)*δ^3*sratio/λ2^7 - 2*exp((-1/2)*δ^2/λ2^2)*δ^5*sratio/λ2^9
    end

    if dorder_view==[2, 0, 1, 1]
        func = 4*exp((-1/2)*δ^2/λ2^2)*sratio/λ2^3 - 10*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio/λ2^5 + 2*exp((-1/2)*δ^2/λ2^2)*δ^4*sratio/λ2^7
    end

    if dorder_view==[1, 0, 1, 1]
        func = 4*exp((-1/2)*δ^2/λ2^2)*δ*sratio/λ2^3 - 2*exp((-1/2)*δ^2/λ2^2)*δ^3*sratio/λ2^5
    end

    if dorder_view==[0, 0, 1, 1]
        func = 2*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio/λ2^3
    end

    if dorder_view==[6, 1, 0, 1]
        func = 0
    end

    if dorder_view==[5, 1, 0, 1]
        func = 0
    end

    if dorder_view==[4, 1, 0, 1]
        func = 0
    end

    if dorder_view==[3, 1, 0, 1]
        func = 0
    end

    if dorder_view==[2, 1, 0, 1]
        func = 0
    end

    if dorder_view==[1, 1, 0, 1]
        func = 0
    end

    if dorder_view==[0, 1, 0, 1]
        func = 0
    end

    if dorder_view==[6, 0, 0, 1]
        func = -30*exp((-1/2)*δ^2/λ2^2)*sratio/λ2^6 + 90*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio/λ2^8 - 30*exp((-1/2)*δ^2/λ2^2)*δ^4*sratio/λ2^10 + 2*exp((-1/2)*δ^2/λ2^2)*δ^6*sratio/λ2^12
    end

    if dorder_view==[5, 0, 0, 1]
        func = -30*exp((-1/2)*δ^2/λ2^2)*δ*sratio/λ2^6 + 20*exp((-1/2)*δ^2/λ2^2)*δ^3*sratio/λ2^8 - 2*exp((-1/2)*δ^2/λ2^2)*δ^5*sratio/λ2^10
    end

    if dorder_view==[4, 0, 0, 1]
        func = 6*exp((-1/2)*δ^2/λ2^2)*sratio/λ2^4 - 12*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio/λ2^6 + 2*exp((-1/2)*δ^2/λ2^2)*δ^4*sratio/λ2^8
    end

    if dorder_view==[3, 0, 0, 1]
        func = 6*exp((-1/2)*δ^2/λ2^2)*δ*sratio/λ2^4 - 2*exp((-1/2)*δ^2/λ2^2)*δ^3*sratio/λ2^6
    end

    if dorder_view==[2, 0, 0, 1]
        func = -2*exp((-1/2)*δ^2/λ2^2)*sratio/λ2^2 + 2*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio/λ2^4
    end

    if dorder_view==[1, 0, 0, 1]
        func = -2*exp((-1/2)*δ^2/λ2^2)*δ*sratio/λ2^2
    end

    if dorder_view==[0, 0, 0, 1]
        func = 2*exp((-1/2)*δ^2/λ2^2)*sratio
    end

    if dorder_view==[6, 0, 2, 0]
        func = -630*exp((-1/2)*δ^2/λ2^2)*sratio^2/λ2^8 + 3465*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio^2/λ2^10 - 2520*exp((-1/2)*δ^2/λ2^2)*δ^4*sratio^2/λ2^12 + 546*exp((-1/2)*δ^2/λ2^2)*δ^6*sratio^2/λ2^14 - 42*exp((-1/2)*δ^2/λ2^2)*δ^8*sratio^2/λ2^16 + exp((-1/2)*δ^2/λ2^2)*δ^10*sratio^2/λ2^18
    end

    if dorder_view==[5, 0, 2, 0]
        func = -630*exp((-1/2)*δ^2/λ2^2)*δ*sratio^2/λ2^8 + 945*exp((-1/2)*δ^2/λ2^2)*δ^3*sratio^2/λ2^10 - 315*exp((-1/2)*δ^2/λ2^2)*δ^5*sratio^2/λ2^12 + 33*exp((-1/2)*δ^2/λ2^2)*δ^7*sratio^2/λ2^14 - exp((-1/2)*δ^2/λ2^2)*δ^9*sratio^2/λ2^16
    end

    if dorder_view==[4, 0, 2, 0]
        func = 60*exp((-1/2)*δ^2/λ2^2)*sratio^2/λ2^6 - 285*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio^2/λ2^8 + 165*exp((-1/2)*δ^2/λ2^2)*δ^4*sratio^2/λ2^10 - 25*exp((-1/2)*δ^2/λ2^2)*δ^6*sratio^2/λ2^12 + exp((-1/2)*δ^2/λ2^2)*δ^8*sratio^2/λ2^14
    end

    if dorder_view==[3, 0, 2, 0]
        func = 60*exp((-1/2)*δ^2/λ2^2)*δ*sratio^2/λ2^6 - 75*exp((-1/2)*δ^2/λ2^2)*δ^3*sratio^2/λ2^8 + 18*exp((-1/2)*δ^2/λ2^2)*δ^5*sratio^2/λ2^10 - exp((-1/2)*δ^2/λ2^2)*δ^7*sratio^2/λ2^12
    end

    if dorder_view==[2, 0, 2, 0]
        func = -6*exp((-1/2)*δ^2/λ2^2)*sratio^2/λ2^4 + 27*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio^2/λ2^6 - 12*exp((-1/2)*δ^2/λ2^2)*δ^4*sratio^2/λ2^8 + exp((-1/2)*δ^2/λ2^2)*δ^6*sratio^2/λ2^10
    end

    if dorder_view==[1, 0, 2, 0]
        func = -6*exp((-1/2)*δ^2/λ2^2)*δ*sratio^2/λ2^4 + 7*exp((-1/2)*δ^2/λ2^2)*δ^3*sratio^2/λ2^6 - exp((-1/2)*δ^2/λ2^2)*δ^5*sratio^2/λ2^8
    end

    if dorder_view==[0, 0, 2, 0]
        func = -3*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio^2/λ2^4 + exp((-1/2)*δ^2/λ2^2)*δ^4*sratio^2/λ2^6
    end

    if dorder_view==[6, 1, 1, 0]
        func = 0
    end

    if dorder_view==[5, 1, 1, 0]
        func = 0
    end

    if dorder_view==[4, 1, 1, 0]
        func = 0
    end

    if dorder_view==[3, 1, 1, 0]
        func = 0
    end

    if dorder_view==[2, 1, 1, 0]
        func = 0
    end

    if dorder_view==[1, 1, 1, 0]
        func = 0
    end

    if dorder_view==[0, 1, 1, 0]
        func = 0
    end

    if dorder_view==[6, 0, 1, 0]
        func = 90*exp((-1/2)*δ^2/λ2^2)*sratio^2/λ2^7 - 375*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio^2/λ2^9 + 195*exp((-1/2)*δ^2/λ2^2)*δ^4*sratio^2/λ2^11 - 27*exp((-1/2)*δ^2/λ2^2)*δ^6*sratio^2/λ2^13 + exp((-1/2)*δ^2/λ2^2)*δ^8*sratio^2/λ2^15
    end

    if dorder_view==[5, 0, 1, 0]
        func = 90*exp((-1/2)*δ^2/λ2^2)*δ*sratio^2/λ2^7 - 95*exp((-1/2)*δ^2/λ2^2)*δ^3*sratio^2/λ2^9 + 20*exp((-1/2)*δ^2/λ2^2)*δ^5*sratio^2/λ2^11 - exp((-1/2)*δ^2/λ2^2)*δ^7*sratio^2/λ2^13
    end

    if dorder_view==[4, 0, 1, 0]
        func = -12*exp((-1/2)*δ^2/λ2^2)*sratio^2/λ2^5 + 39*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio^2/λ2^7 - 14*exp((-1/2)*δ^2/λ2^2)*δ^4*sratio^2/λ2^9 + exp((-1/2)*δ^2/λ2^2)*δ^6*sratio^2/λ2^11
    end

    if dorder_view==[3, 0, 1, 0]
        func = -12*exp((-1/2)*δ^2/λ2^2)*δ*sratio^2/λ2^5 + 9*exp((-1/2)*δ^2/λ2^2)*δ^3*sratio^2/λ2^7 - exp((-1/2)*δ^2/λ2^2)*δ^5*sratio^2/λ2^9
    end

    if dorder_view==[2, 0, 1, 0]
        func = 2*exp((-1/2)*δ^2/λ2^2)*sratio^2/λ2^3 - 5*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio^2/λ2^5 + exp((-1/2)*δ^2/λ2^2)*δ^4*sratio^2/λ2^7
    end

    if dorder_view==[1, 0, 1, 0]
        func = 2*exp((-1/2)*δ^2/λ2^2)*δ*sratio^2/λ2^3 - exp((-1/2)*δ^2/λ2^2)*δ^3*sratio^2/λ2^5
    end

    if dorder_view==[0, 0, 1, 0]
        func = exp((-1/2)*δ^2/λ2^2)*δ^2*sratio^2/λ2^3
    end

    if dorder_view==[6, 2, 0, 0]
        func = -630*exp((-1/2)*δ^2/λ1^2)/λ1^8 + 3465*exp((-1/2)*δ^2/λ1^2)*δ^2/λ1^10 - 2520*exp((-1/2)*δ^2/λ1^2)*δ^4/λ1^12 + 546*exp((-1/2)*δ^2/λ1^2)*δ^6/λ1^14 - 42*exp((-1/2)*δ^2/λ1^2)*δ^8/λ1^16 + exp((-1/2)*δ^2/λ1^2)*δ^10/λ1^18
    end

    if dorder_view==[5, 2, 0, 0]
        func = -630*exp((-1/2)*δ^2/λ1^2)*δ/λ1^8 + 945*exp((-1/2)*δ^2/λ1^2)*δ^3/λ1^10 - 315*exp((-1/2)*δ^2/λ1^2)*δ^5/λ1^12 + 33*exp((-1/2)*δ^2/λ1^2)*δ^7/λ1^14 - exp((-1/2)*δ^2/λ1^2)*δ^9/λ1^16
    end

    if dorder_view==[4, 2, 0, 0]
        func = 60*exp((-1/2)*δ^2/λ1^2)/λ1^6 - 285*exp((-1/2)*δ^2/λ1^2)*δ^2/λ1^8 + 165*exp((-1/2)*δ^2/λ1^2)*δ^4/λ1^10 - 25*exp((-1/2)*δ^2/λ1^2)*δ^6/λ1^12 + exp((-1/2)*δ^2/λ1^2)*δ^8/λ1^14
    end

    if dorder_view==[3, 2, 0, 0]
        func = 60*exp((-1/2)*δ^2/λ1^2)*δ/λ1^6 - 75*exp((-1/2)*δ^2/λ1^2)*δ^3/λ1^8 + 18*exp((-1/2)*δ^2/λ1^2)*δ^5/λ1^10 - exp((-1/2)*δ^2/λ1^2)*δ^7/λ1^12
    end

    if dorder_view==[2, 2, 0, 0]
        func = -6*exp((-1/2)*δ^2/λ1^2)/λ1^4 + 27*exp((-1/2)*δ^2/λ1^2)*δ^2/λ1^6 - 12*exp((-1/2)*δ^2/λ1^2)*δ^4/λ1^8 + exp((-1/2)*δ^2/λ1^2)*δ^6/λ1^10
    end

    if dorder_view==[1, 2, 0, 0]
        func = -6*exp((-1/2)*δ^2/λ1^2)*δ/λ1^4 + 7*exp((-1/2)*δ^2/λ1^2)*δ^3/λ1^6 - exp((-1/2)*δ^2/λ1^2)*δ^5/λ1^8
    end

    if dorder_view==[0, 2, 0, 0]
        func = -3*exp((-1/2)*δ^2/λ1^2)*δ^2/λ1^4 + exp((-1/2)*δ^2/λ1^2)*δ^4/λ1^6
    end

    if dorder_view==[6, 1, 0, 0]
        func = 90*exp((-1/2)*δ^2/λ1^2)/λ1^7 - 375*exp((-1/2)*δ^2/λ1^2)*δ^2/λ1^9 + 195*exp((-1/2)*δ^2/λ1^2)*δ^4/λ1^11 - 27*exp((-1/2)*δ^2/λ1^2)*δ^6/λ1^13 + exp((-1/2)*δ^2/λ1^2)*δ^8/λ1^15
    end

    if dorder_view==[5, 1, 0, 0]
        func = 90*exp((-1/2)*δ^2/λ1^2)*δ/λ1^7 - 95*exp((-1/2)*δ^2/λ1^2)*δ^3/λ1^9 + 20*exp((-1/2)*δ^2/λ1^2)*δ^5/λ1^11 - exp((-1/2)*δ^2/λ1^2)*δ^7/λ1^13
    end

    if dorder_view==[4, 1, 0, 0]
        func = -12*exp((-1/2)*δ^2/λ1^2)/λ1^5 + 39*exp((-1/2)*δ^2/λ1^2)*δ^2/λ1^7 - 14*exp((-1/2)*δ^2/λ1^2)*δ^4/λ1^9 + exp((-1/2)*δ^2/λ1^2)*δ^6/λ1^11
    end

    if dorder_view==[3, 1, 0, 0]
        func = -12*exp((-1/2)*δ^2/λ1^2)*δ/λ1^5 + 9*exp((-1/2)*δ^2/λ1^2)*δ^3/λ1^7 - exp((-1/2)*δ^2/λ1^2)*δ^5/λ1^9
    end

    if dorder_view==[2, 1, 0, 0]
        func = 2*exp((-1/2)*δ^2/λ1^2)/λ1^3 - 5*exp((-1/2)*δ^2/λ1^2)*δ^2/λ1^5 + exp((-1/2)*δ^2/λ1^2)*δ^4/λ1^7
    end

    if dorder_view==[1, 1, 0, 0]
        func = 2*exp((-1/2)*δ^2/λ1^2)*δ/λ1^3 - exp((-1/2)*δ^2/λ1^2)*δ^3/λ1^5
    end

    if dorder_view==[0, 1, 0, 0]
        func = exp((-1/2)*δ^2/λ1^2)*δ^2/λ1^3
    end

    if dorder_view==[6, 0, 0, 0]
        func = -15*exp((-1/2)*δ^2/λ1^2)/λ1^6 + 45*exp((-1/2)*δ^2/λ1^2)*δ^2/λ1^8 - 15*exp((-1/2)*δ^2/λ1^2)*δ^4/λ1^10 + exp((-1/2)*δ^2/λ1^2)*δ^6/λ1^12 - 15*exp((-1/2)*δ^2/λ2^2)*sratio^2/λ2^6 + 45*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio^2/λ2^8 - 15*exp((-1/2)*δ^2/λ2^2)*δ^4*sratio^2/λ2^10 + exp((-1/2)*δ^2/λ2^2)*δ^6*sratio^2/λ2^12
    end

    if dorder_view==[5, 0, 0, 0]
        func = -15*exp((-1/2)*δ^2/λ1^2)*δ/λ1^6 + 10*exp((-1/2)*δ^2/λ1^2)*δ^3/λ1^8 - exp((-1/2)*δ^2/λ1^2)*δ^5/λ1^10 - 15*exp((-1/2)*δ^2/λ2^2)*δ*sratio^2/λ2^6 + 10*exp((-1/2)*δ^2/λ2^2)*δ^3*sratio^2/λ2^8 - exp((-1/2)*δ^2/λ2^2)*δ^5*sratio^2/λ2^10
    end

    if dorder_view==[4, 0, 0, 0]
        func = 3*exp((-1/2)*δ^2/λ1^2)/λ1^4 - 6*exp((-1/2)*δ^2/λ1^2)*δ^2/λ1^6 + exp((-1/2)*δ^2/λ1^2)*δ^4/λ1^8 + 3*exp((-1/2)*δ^2/λ2^2)*sratio^2/λ2^4 - 6*exp((-1/2)*δ^2/λ2^2)*δ^2*sratio^2/λ2^6 + exp((-1/2)*δ^2/λ2^2)*δ^4*sratio^2/λ2^8
    end

    if dorder_view==[3, 0, 0, 0]
        func = 3*exp((-1/2)*δ^2/λ1^2)*δ/λ1^4 - exp((-1/2)*δ^2/λ1^2)*δ^3/λ1^6 + 3*exp((-1/2)*δ^2/λ2^2)*δ*sratio^2/λ2^4 - exp((-1/2)*δ^2/λ2^2)*δ^3*sratio^2/λ2^6
    end

    if dorder_view==[2, 0, 0, 0]
        func = -exp((-1/2)*δ^2/λ1^2)/λ1^2 + exp((-1/2)*δ^2/λ1^2)*δ^2/λ1^4 - exp((-1/2)*δ^2/λ2^2)*sratio^2/λ2^2 + exp((-1/2)*δ^2/λ2^2)*δ^2*sratio^2/λ2^4
    end

    if dorder_view==[1, 0, 0, 0]
        func = -exp((-1/2)*δ^2/λ1^2)*δ/λ1^2 - exp((-1/2)*δ^2/λ2^2)*δ*sratio^2/λ2^2
    end

    if dorder_view==[0, 0, 0, 0]
        func = exp((-1/2)*δ^2/λ2^2)*sratio^2 + exp((-1/2)*δ^2/λ1^2)
    end

    dorder[2] = dorder2  # resetting dorder[2]
    return powers_of_negative_one(dorder2) * float(func)  # correcting for amount of t2 derivatives

end

return se_se_kernel, 3  # the function handle and the number of kernel hyperparameters

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
export Melody

using ..NoteSequences: DEFAULT_STEPS_PER_BAR, DEFAULT_STEPS_PER_QUARTER
using ..NoteSequences: DEFAULT_QPM, DEFAULT_TPQ
using MacroTools: @forward

"""
    Melody <: Any

Melody is an intermediate representation for monophonic musical sequences.
Melody events are integers in the range [-2, 127].
The negative values are `MELODY_NO_EVENt`, `MELODY_NOTE_OFF` and they are special events.
Note on events are the non-negative values from [0, 127] and they represent a midi pitch.

A melody note starts at a midi pitch value (non-negative). The note is sustained through
the following `MELODY_NO_EVENT`s until another midi pitch is reached or a `MELODY_NOTE_OFF`
event is reached.

## Fields
* `events::Vector{Int}`: The monophonic melody events in the sequence.
* `steps_per_bar::Int`: Number of steps per bar of music.
* `steps_per_quarter::Int`: Number of steps per quarter note.
* `startstep::Int`: The first step of the melody in the sequence.
* `endstep::Int`: The last step of the melody in the sequence.
"""
mutable struct Melody
    events::Vector{Int}
    steps_per_bar::Int
    steps_per_quarter::Int
    startstep::Int
    endstep::Int

    function Melody(;events::Vector{Int}=Int[],
                     startstep::Int=0,
                     steps_per_bar::Int=DEFAULT_STEPS_PER_BAR,
                     steps_per_quarter::Int=DEFAULT_STEPS_PER_QUARTER)

        if !isempty(events)
            for event in events
                isoutofrange(event) &&
                    throw(error("Melody event is out of range"))
            end

            # Up to the first event, replace MELODY_NOTE_OFF with MELODY_NO_EVENT
            # to prevent unwanted note-offs.
            for (i, event) in enumerate(events)
                event ∉ [MELODY_NO_EVENT, MELODY_NOTE_OFF] && break
                events[i] = MELODY_NO_EVENT
            end
        end

        endstep = startstep + length(events)
        new(events, steps_per_bar, steps_per_quarter, startstep, endstep)
    end
end

@forward Melody.events Base.getindex, Base.length, Base.iterate, Base.lastindex

function Base.push!(melody::Melody, event::Int)
    isoutofrange(event) && throw(error("Melody event is out of range"))
    push!(melody.events, event)
    melody.endstep += 1
    melody
end

function Base.append!(melody::Melody, events::Vector{Int})
    append!(melody.events, events)
    melody.endstep += length(events)
    melody
end

Base.setindex!(m::Melody, event::Int, idx::Int) = setindex!(m.events, event, idx)

"""
    Melody(qns::NoteSequence,
                search_start_step::Int=0,
                gapbars::Int=1,
                ignore_polyphony::Bool=false,
                pad_end::Bool=false,
                instrument::Int=1)

Convert the given quantized notesequence to a melody. The melody is extracted starting from
`search_start_step`. The melody extraction will be ended when there is a silence for more than
`gapbars` amount of bars in music.

If `ignore_polyphony` is true, only the highest pitch is kept when multiple pitches are encountered.
Otherwise, an error is thrown.

0 velocity notes are ignored when extracting the melody. If `pad_end` is true, the melody is padded with
`MELODY_NO_EVENt` until the end of the bar.
"""
function Melody(qns::NoteSequence,
                search_start_step::Int=0,
                gapbars::Int=1,
                ignore_polyphony::Bool=false,
                pad_end::Bool=false,
                instrument::Int=1)

    !qns.isquantized && throw(error("NoteSequence must be relative quantized."))

    quarters_per_beat = 4 / qns.timesignatures[1].denominator
    quarters_per_bar = quarters_per_beat * qns.timesignatures[1].numerator
    steps_per_bar = quarters_per_bar * qns.steps_per_quarter
    steps_per_bar % 1 != 0 && throw(error("Steps per bar must be an integer"))

    melody = Melody()
    melody.steps_per_bar = round(steps_per_bar)
    melody.steps_per_quarter = qns.steps_per_quarter

    notes = [note for note in qns.notes
             if note.instrument == instrument && note.start_time >= search_start_step]

    isempty(notes) && return

    sort!(notes, by=note -> (note.start_time, -note.pitch))

    initial_distance = (notes[1].start_time - search_start_step)
    melody_start_step = notes[1].start_time - initial_distance % melody.steps_per_bar

    for note in notes
        note.velocity == 0 && continue

        note_start_step = note.start_time - melody_start_step
        note_end_step = note.end_time - melody_start_step

        if isempty(melody.events)
            addnote!(melody, note.pitch, note_start_step, note_end_step)
            continue
        end

        laston, lastoff = last_onoff_events(melody)

        on_distance = note_start_step - laston
        off_distance = note_start_step - lastoff

        if on_distance == 0
            if ignore_polyphony
                continue # Keeps the highest note
            else
                throw(error("More than one note is found at the same time."))
            end
        elseif on_distance < 0
            throw(error("Notes not in ascending order. This is caused due to polyphonic melody."))
        end

        gapsteps = gapbars * melody.steps_per_bar

        # End the melody if there is a silence of `gap_bars` or More
        !isempty(melody) && off_distance >= gapsteps && break

        addnote!(melody, note.pitch, note_start_step, note_end_step)
    end

    isempty(melody) && return

    melody.startstep = melody_start_step
    melody[end] == MELODY_NOTE_OFF && pop!(melody.events)

    melodylength = length(melody)
    if pad_end
        melodylength += (-melodylength % melody.steps_per_bar)
    end
    setlength!(melody, melodylength)

    melody
end

"""
    getnotesequence(melody::Melody;
                    velocity::Int=100,
                    instrument::Int=1,
                    program::Int=0,
                    sequence_start_time::Int=0,
                    qpm::Float64=DEFAULT_QPM)

Convert the given melody to a `NoteSequence`. The notes events will have the same `velocity`,
`program` and `instrument`. The melody will start at `sequence_start_time` in the sequence
and will have a tempo of `qpm` (quarter notes per minute).
"""
function getnotesequence(melody::Melody;
                         velocity::Int=100,
                         instrument::Int=1,
                         program::Int=0,
                         sequence_start_time::Int=0,
                         qpm::Float64=DEFAULT_QPM)

    seconds_per_step = 60 / qpm / melody.steps_per_quarter
    ns = NoteSequence()
    push!(ns.tempos, NoteSequences.Tempo(0, qpm))
    ns.tpq = DEFAULT_TPQ

    sequence_start_time += melody.startstep * seconds_per_step
    note_is_playing = false
    local current_seqnote

    for (step, note) in enumerate(melody)
        seconds = step * seconds_per_step + sequence_start_time
        new_note_start_time = second2tick(seconds)

        if MIN_MIDI_PITCH <= note <= MAX_MIDI_PITCH
            if note_is_playing
                # End the sustained note
                current_seqnote.end_time = new_note_start_time
            end

            current_seqnote = NoteSequences.SeqNote(note, velocity,
                                                    new_note_start_time, 0,
                                                    instrument, program)

            push!(ns.notes, current_seqnote)
            note_is_playing = true

        elseif note == MELODY_NOTE_OFF
            if note_is_playing
                current_seqnote.end_time = new_note_start_time
                note_is_playing = false
            end
        end
    end

    if note_is_playing
        current_seqnote.end_time = second2tick(length(melody) * seconds_per_step + sequence_start_time)
    end

    if !isempty(ns.notes)
        ns.total_time = ns.notes[end].end_time
    end

    ns
end

"""
    setlength!(melody::Melody, steps::Int)

Set the length of the melody to `steps`. If the given `steps` is greater than
the length of the melody, the melody will be padded with `MELODY_NO_EVENT`s to the right.
If `steps` is smaller than the length of the melody, the melody will be truncated.
"""
function setlength!(melody::Melody, steps::Int)
    oldlength = length(melody)

    if steps > length(melody)
        padding = fill(MELODY_NO_EVENT, (steps - length(melody)))
        append!(melody, padding)
    else
        melody.events = melody.events[1:steps]
    end

    melody.endstep = melody.startstep + steps

    if steps > oldlength
        for i = oldlength:-1:1
            if melody[i] == MELODY_NOTE_OFF
                break
            elseif melody[i] != MELODY_NO_EVENT
                melody[oldlength] = MELODY_NOTE_OFF
            end
        end
    end

    melody
end

"""
    addnote!(melody::Melody, pitch::Int, startstep::Int, endstep::Int)

Adds the given note to the melody at `startstep` and `endstep`.
The `endstep` will be set to a `MELODY_NOTE_OFF`.
"""
function addnote!(melody::Melody, pitch::Int, startstep::Int, endstep::Int)
    startstep >= endstep && throw(error("End step must be greater than start step."))
    setlength!(melody, endstep)
    melody[startstep] = pitch
    melody[endstep] = MELODY_NOTE_OFF

    for i = (startstep + 1):endstep
        melody[i] = MELODY_NO_EVENT
    end

    melody
end

"""
    last_onoff_events(melody::Melody)

Return the indices of the last midi pitch and `MELODY_NOTE_OFF` in the sequence.
"""
function last_onoff_events(melody::Melody)
    lastoff = length(melody)
    for i = (length(melody)-1):-1:-1
        if melody[i] == MELODY_NOTE_OFF
            lastoff = i
        end
        melody[i] >= MIN_MIDI_PITCH && return (i, lastoff)
    end
end

isoutofrange(event::Int) = !(MIN_MELODY_EVENT <= event <= MAX_MELODY_EVENT)


text:
abstract type PeriodicOrbitAlgorithm end
"""
	guessFromHopf(br, ind_hopf, eigsolver::AbstractEigenSolver, M, amplitude; phase = 0)

This function returns several useful quantities regarding a Hopf bifurcation point. More precisely, it returns:
- the parameter value at which a Hopf bifurcation occurs
- the period of the bifurcated periodic orbit
- a guess for the bifurcated periodic orbit
- the equilibrium at the Hopf bifurcation point
- the eigenvector at the Hopf bifurcation point.

The arguments are
- `br`: the continuation branch which lists the Hopf bifurcation points
- `ind_hopf`: index of the bifurcation branch, as in `br.specialpoint`
- `eigsolver`: the eigen solver used to find the eigenvectors
- `M` number of time slices in the periodic orbit guess
- `amplitude`: amplitude of the periodic orbit guess
"""
function guessFromHopf(br, ind_hopf, eigsolver::AbstractEigenSolver, M, amplitude; phase = 0)
	hopfpoint = HopfPoint(br, ind_hopf)
	specialpoint = br.specialpoint[ind_hopf]

	# parameter value at the Hopf point
	p_hopf = hopfpoint.p[1]

	# frequency at the Hopf point
	ωH  = hopfpoint.p[end] |> abs

	# vec_hopf is the eigenvector for the eigenvalues iω
	vec_hopf = geteigenvector(eigsolver, br.eig[specialpoint.idx][2], specialpoint.ind_ev-1)
	vec_hopf ./=  norm(vec_hopf)

	 orbitguess = [real.(hopfpoint.u .+ amplitude .* vec_hopf .* exp(-2pi * complex(0, 1) .* (ii/(M-1) - phase))) for ii=0:M-1]

	return p_hopf, 2pi/ωH, orbitguess, hopfpoint, vec_hopf
end
####################################################################################################
# Amplitude of the u component of the cycle
amplitude(x::AbstractMatrix, n) =  maximum(x[1:n, :]) - minimum(x[1:n, :])

function amplitude(x::AbstractVector, n, M; ratio = 1)
	xc = reshape(x[1:end-1], ratio * n, M)
	amplitude(xc, n)
end

function maximumPOTrap(x::AbstractVector, n, M; ratio = 1)
	xc = reshape(x[1:end-1], ratio * n, M)
	maximum(x[1:n, :])
end


text:
struct ClusterTree{T}
    npar::Int64
    parindices::Vector{Int64}
    root::Cluster{T}
end

function ClusterTree(particles::Particles; n, threshold, stretch=SVector(1.0,1.0,1.0)) where {T}
    npar = particles.npar
    parindices = [1:npar;]
    root = make_cluster(particles, parindices, 1, npar, 1; n=n, threshold=threshold, stretch=stretch)
    return ClusterTree(npar, parindices, root)
end

function print_bbox(ct::ClusterTree{T}) where {T}
    print_bbox(ct.root)
end


text:
function x2c(ψ::Array{Complex{Float64},1},tinfo::ProjectedGPE.Tinfo)
    return tinfo.Tx'*(tinfo.W.*ψ)
end

function x2c(ψ::Array{Complex{Float64},2},tinfo::ProjectedGPE.Tinfo)
    return tinfo.Tx'*(tinfo.W.*ψ)*tinfo.Ty
end

function x2c(ψ::Array{Complex{Float64},3},tinfo::ProjectedGPE.Tinfo)
    sa = size(ψ);sx = size(tinfo.Tx');sy = size(tinfo.Ty');sz = size(tinfo.Tz')
    Ax = reshape(tinfo.W.*ψ,(sa[1],sa[2]*sa[3]))
    Ax = reshape(tinfo.Tx'*Ax,(sx[1],sa[2],sa[3]))
    Ax = permutedims(Ax,(2,3,1));sa = size(Ax)

    Ay = reshape(Ax,(sa[1],sa[2]*sa[3]))
    Ay = reshape(tinfo.Ty'*Ay,(sy[1],sa[2],sa[3]))
    Ay = permutedims(Ay,(2,3,1));sa = size(Ay)

    Az = reshape(Ay,(sa[1],sa[2]*sa[3]))
    Az = reshape(tinfo.Tz'*Az,(sz[1],sa[2],sa[3]))
    return permutedims(Az,(2,3,1))

end


text:
using Mitosis, LinearAlgebra, Test, Statistics
m = [1.0, 0.5]
K = 10000
Q = Matrix(1.0I, 2, 2)
G = Gaussian(μ=m, Σ=Q)

@testset "Gauss" begin
    @test G === Gaussian{keys(G)}(m, Q)

    @test m === mean(G)
    @test Q === cov(G)

    @test 10/sqrt(K) > norm(m - mean(rand(G) for x in 1:K))
    @test 10/sqrt(K) > norm(Q - cov([rand(G) for x in 1:K]))
end

p = WGaussian(μ=m, Σ=Q, c=0.0)
p2 = convert(WGaussian{(:μ,:Σ,:c)}, convert(WGaussian{(:F,:Γ,:c)}, p))
@test p ≈ p2

p = Gaussian(μ=m, Σ=Q)
p2 = convert(Gaussian{(:μ,:Σ)}, convert(Gaussian{(:F,:Γ)}, p))
@test p ≈ p2


text:
#
using DiffMesh, LinearAlgebra

n = 100

dl = rand(n-1)
dd = rand(n)
du = rand(n-1)

Td = Tridiagonal(dl,dd,du)
Tt = Tdg(dl,dd,du)

u = rand(n)

@test Td *u ≈ Tt *u atol=1e-12
@test Td'*u ≈ Tt'*u atol=1e-12


text:
module test_julia_core

using Test

module A
_val = nothing
end

@test A._val === nothing

Core.eval(A, :(_val = 1))

@test A._val == 1

end # module test_julia_core

module test_julia_core_compiler

using Test
using Core.Compiler: ⊑

@test Int ⊑ Integer

end # module test_julia_core_compiler


text:
# Practical Machine learning
# Bayesian learning - Naive Bayes example
# Chapter 9

import StatsBase.fit
import StatsBase.predict
using Distributions

#### core naive Bayes functions ####

function ensure_data_size(X, y)
    @assert(size(X, 2) == length(y),
            "Number of observations in X ($(size(X, 2))) is not equal to " *
            "number of class labels in y ($(length(y)))")
end

function logprob_c{C}(m::NBModel, c::C)
    return m.c_counts[c] / m.n_obs
end

# predict log probabilities for all classes
function predict_logprobs{V<:Number}(m::NBModel, x::Vector{V})
    C = eltype(m.c_counts)[1]
    logprobs = Dict{C, Float64}()
    for c in keys(m.c_counts)
        logprobs[c] = logprob_c(m, c) + logprob_x_given_c(m, x, c)
    end
    return keys(logprobs), values(logprobs)
end

# predict log probabilities for all classes
function predict_logprobs{V<:Number}(m::NBModel, X::Matrix{V})
    C = eltype(m.c_counts)[1]
    logprobs_per_class = Dict{C, Vector{Float64}}()
    for c in keys(m.c_counts)
        logprobs_per_class[c] = logprob_c(m, c) + logprob_x_given_c(m, X, c)
    end
    return (collect(keys(logprobs_per_class)),
            hcat(collect(values(logprobs_per_class))...)')
end

# preditct logprobs, return tuples of predicted class and its logprob
function predict_proba{V<:Number}(m::NBModel, X::Matrix{V})
    C = eltype(m.c_counts)[1]
    classes, logprobs = predict_logprobs(m, X)
    predictions = Array((C, Float64), size(X, 2))
    for j=1:size(X, 2)
        maxprob_idx = indmax(logprobs[:, j])
        c = classes[maxprob_idx]
        logprob = logprobs[maxprob_idx, j]
        predictions[j] = (c, logprob)
    end
    return predictions
end

function predict{V<:Number}(m::NBModel, X::Matrix{V})
    return [k for (k,v) in predict_proba(m, X)]
end

#####  Multinomial Naive Bayes  #####

function fit{C}(m::MultinomialNB, X::Matrix{Int64}, y::Vector{C})
    ensure_data_size(X, y)
    for j=1:size(X, 2)
        c = y[j]
        m.c_counts[c] += 1
        m.x_counts[c] .+= X[:, j]
        m.x_totals += X[:, j]
        m.n_obs += 1
    end
    return m
end

# Calculate log P(x|C)
function logprob_x_given_c{C}(m::MultinomialNB, x::Vector{Int64}, c::C)
    x_priors_for_c = m.x_counts[c] ./ m.x_totals
    x_probs_given_c = x_priors_for_c .^ x
    logprob = sum(log(x_probs_given_c))
    return logprob
end

# Calculate log P(x|C)
function logprob_x_given_c{C}(m::MultinomialNB, X::Matrix{Int64}, c::C)
    x_priors_for_c = m.x_counts[c] ./ m.x_totals
    x_probs_given_c = x_priors_for_c .^ X
    logprob = sum(log(x_probs_given_c), 1)
    return squeeze(logprob, 1)
end

######  Gaussian Naive Bayes  #######

function fit{C}(m::GaussianNB, X::Matrix{Float64}, y::Vector{C})
    ensure_data_size(X, y)
    # updatestats(m.dstats, X)
    # m.gaussian = MvNormal(mean(m.dstats), cov(m.dstats))
    # m.n_obs = m.dstats.n_obs
    n_vars = size(X, 1)
    for j=1:size(X, 2)
        c = y[j]
        m.c_counts[c] += 1
        updatestats(m.c_stats[c], reshape(X[:, j], n_vars, 1))
        # m.x_counts[c] .+= X[:, j]
        # m.x_totals += X[:, j]
        m.n_obs += 1
    end
    # precompute distributions for each class
    for c in keys(m.c_counts)
        m.gaussians[c] = MvNormal(mean(m.c_stats[c]), cov(m.c_stats[c]))
    end
    return m
end

# Calculate log P(x|C)
function logprob_x_given_c{C}(m::GaussianNB, x::Vector{Float64}, c::C)
    return logpdf(m.gaussians[c], x)
end

# Calculate log P(x|C)
function logprob_x_given_c{C}(m::GaussianNB, X::Matrix{Float64}, c::C)
    ## x_priors_for_c = m.x_counts[c] ./ m.x_totals
    ## x_probs_given_c = x_priors_for_c .^ x
    ## logprob = sum(log(x_probs_given_c))
    ## return logprob
    return logpdf(m.gaussians[c], X)
end


text:
l = @latexify dummyfunc(x; y=1, z=3) = x^2/y + z
@test l == raw"$\mathrm{dummyfunc}\left( x; y = 1, z = 3 \right) = \frac{x^{2}}{y} + z$"

@test_throws UndefVarError dummyfunc(1.)

l2 = @latexrun dummyfunc2(x; y=1, z=3) = x^2/y + z
@test l2 == raw"$\mathrm{dummyfunc2}\left( x; y = 1, z = 3 \right) = \frac{x^{2}}{y} + z$"

@test dummyfunc2(1.) == 4

l3 = @latexify dummyfunc2(x::Number; y=1, z=3) = x^2/y + z
@test l3 == raw"$\mathrm{dummyfunc2}\left( x::Number; y = 1, z = 3 \right) = \frac{x^{2}}{y} + z$"

l4 = @latexify dummyfunc2(::Number; y=1, z=3) = x^2/y + z
@test l4 == raw"$\mathrm{dummyfunc2}\left( ::Number; y = 1, z = 3 \right) = \frac{x^{2}}{y} + z$"

@test latexify(:(@hi(x / y))) == replace(
raw"$\mathrm{@hi}\left( \frac{x}{y} \right)$", "\r\n"=>"\n")


text:
# power analysis functions for t-test

"""
    powerTTest(d::Real = 0, n::Real = 0, alpha::Float64 = 0.05, alternative = "two.sided")

Compute power of a sample of `n` observations to test the effect size `d`
at type I error = `alpha` (default: 0.05).
The effect size `d` is the difference in means divided by the pooled standard deviation. The option
`alternative` is `two.sided` for H₀: μ₁ = μ₂ vs H₁: μ₁ ≠ μ₂, `less` for H₁: μ₁ < μ₂, and `greater`
for H₁: μ₁ > μ₂. `sampletype` is `onesample` for one-sample ttest, `twosample` for two-sample ttest,
and `paired` for two-sample paired ttest.
"""
function powerTTest(;
    d::Real = 0.0,
    n::Real = 0,
    alpha::Float64 = 0.05,
    sampletype::String = "onesample",
    alternative::String = "two")

    check_args(d=d,n=n,alpha=alpha)

    if lowercase(sampletype) in ("onesample","one.sample","one-sample","one sample","paired")
        tsample = 1
    elseif lowercase(sampletype) in ("twosample","two.sample","two-sample","two sample")
        tsample = 2
    end

    if lowercase(alternative) == "less"
        tside = ttside = 1
    elseif lowercase(alternative) in ("two","two.sided","two-sided","two sided")
        tside = ttside = 2
        d = abs(d)
    elseif lowercase(alternative) == "greater"
        ttside = 3
        tside = 1
    end

    ν = ceil(Int64,(n-1)*tsample)
    λ = sqrt(n/tsample)*d
    if ttside == 1
        return cdf(NoncentralT(ν,λ),quantile(TDist(ν),alpha/tside))
    elseif ttside == 2
        qu = cquantile(TDist(ν),alpha/tside)
        return ccdf(NoncentralT(ν,λ),qu) + cdf(NoncentralT(ν,λ),-qu)
    elseif ttside == 3
        return ccdf(NoncentralT(ν,λ),cquantile(TDist(ν),alpha/tside))
    end
end

"""
    samplesizeTTest(d::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the sample size to test the effect size `d` with power > `power` (default: 0.8)
at type I error = `alpha` (default: 0.05).
The effect size `d` is the difference in means divided by the pooled standard deviation. The option
`alternative` is `two.sided` for H₀: μ₁ = μ₂ vs H₁: μ₁ ≠ μ₂ (default), `less` for H₁: μ₁ < μ₂, and `greater`
for H₁: μ₁ > μ₂. `sampletype` is `onesample` for one-sample ttest (default), `twosample` for two-sample ttest,
and `paired` for two-sample paired ttest.
"""
function samplesizeTTest(;
    d::Real = 0.0,
    alpha::Float64 = 0.05,
    power::Float64 = 0.8,
    sampletype::String = "onesample",
    alternative::String = "two")

    check_args(d=d,alpha=alpha,power=power)

    return ceil(Int64,fzero(x->powerTTest(n = x, d = d, alpha = alpha, sampletype = sampletype, alternative = alternative) - power, 2.0, 10.0^7))
end

"""
    effectsizeTTest(n::Real = 0, alpha::Float64 = 0.05, power::Float64 = 0.8, alternative = "two.sided")

Compute the effect size that a sample of `n` can detect with power > `power` (default: 0.8)
at type I error = `alpha` (default: 0.05).
The effect size `d` is the difference in means divided by the pooled standard deviation. The option
`alternative` is `two.sided` for H₀: μ₁ = μ₂ vs H₁: μ₁ ≠ μ₂ (default), `less` for H₁: μ₁ < μ₂, and `greater`
for H₁: μ₁ > μ₂. `sampletype` is `onesample` for one-sample ttest (default), `twosample` for two-sample ttest,
and `paired` for two-sample paired ttest.
"""
function effectsizeTTest(;
    n::Real = 0,
    alpha::Float64 = 0.05,
    power::Float64 = 0.8,
    sampletype::String = "onesample",
    alternative::String = "two"
    )

    check_args(n=n,alpha=alpha,power=power)

    return fzero(x -> powerTTest(n = n, d = x, alpha = alpha, sampletype = sampletype, alternative = alternative) - power,.001,100)
end

"""
    alphaTTest(d::Real = 0, n::Real = 0, power::Float64 = 0.8, alternative = "two.sided")

Compute the probability of type I error that a sample of `n` has in detecting the
effect size `d` with power > `power` (default: 0.8) at type I error = `alpha` (default: 0.05).
The effect size `d` is the difference in means divided by the pooled standard deviation. The option
`alternative` is `two.sided` for H₀: μ₁ = μ₂ vs H₁: μ₁ ≠ μ₂ (default), `less` for H₁: μ₁ < μ₂, and `greater`
for H₁: μ₁ > μ₂. `sampletype` is `onesample` for one-sample ttest (default), `twosample` for two-sample ttest,
and `paired` for two-sample paired ttest.
"""
function alphaTTest(;
    n::Real = 0,
    d::Real = 0.0,
    power::Float64 = 0.8,
    sampletype::String = "onesample",
    alternative::String = "two"
    )

    check_args(d=d,n=n,power=power)

    return fzero(x->powerTTest(n = n, d = d, alpha = x, sampletype = sampletype, alternative = alternative) - power, 1e-10, 1 - 1e-10)
end

"""
    pwr.TTest(d::Real = 0, n::Real = 0, alpha::Float64 = 0.0, power::Float64 = 0.0, alternative = "two.sided")

Compute one of the test parameters such as sample size (`n`),
effect size (`d`), type I error (`alpha`), and power (`power`).
The parameter to be estimated must be set to zero (default).
The effect size `d` is the difference in means divided by the pooled standard deviation. The option
`alternative` is `two.sided` for H₀: μ₁ = μ₂ vs H₁: μ₁ ≠ μ₂ (default), `less` for H₁: μ₁ < μ₂, and `greater`
for H₁: μ₁ > μ₂. `sampletype` is `onesample` for one-sample ttest (default), `twosample` for two-sample ttest,
and `paired` for two-sample paired ttest.
"""
function TTest(;
    n::Real = 0,
    d::Real = 0.0,
    alpha::Float64 = 0.0,
    power::Float64 = 0.0,
    sampletype::String = "onesample",
    alternative::String = "two")

    if sum([x == 0 for x in (n,d,alpha,power)]) != 1
        error("exactly one of n, d, power, and alpha must be zero")
    end

    if power == 0.0
        power = powerTTest(n = n, d = d, alpha = alpha, sampletype = sampletype, alternative = alternative)
    elseif alpha == 0.0
        alpha = alphaTTest(n = n, d = d, power = power, sampletype = sampletype, alternative = alternative)
    elseif d == 0.0
        d = effectsizeTTest(n = n, alpha = alpha, power = power, sampletype = sampletype, alternative = alternative)
    elseif n == 0
        n = samplesizeTTest(d = d, alpha = alpha, power = power, sampletype = sampletype, alternative = alternative)
    end

    stype = Dict("onesample" => "One-sample","one.sample" => "One-sample","one sample" => "One-sample",
        "one-sample" => "One-sample","twosample" => "Two-sample","two.sample" => "Two-sample",
        "two sample" => "Two-sample","two-sample" => "Two-sample","paired" => "Paired")
    alt = Dict("two" => "two-sided", "two.sided" => "two-sided", "less" => "less", "greater" => "greater")

    note = ""
    if alternative == "paired"
        note = "`n` is number of pairs"
    elseif alternative in ("two","two.sided","two-sided","two sided")
        note = "`n` is number in each group"
    end

    return htest(
        string(stype[sampletype]," t-test power calculation"),
        OrderedDict(
            "n" => n,
            "d" => d,
            "alpha" => alpha,
            "power" => power,
            "sampletype" => stype[sampletype],
            "alternative" => alt[alternative],
            "note" => note)
        )
end

max_stars_repo_name: mwsohn/pwr

max_issues_repo_name: mwsohn/pwr

max_forks_repo_name: mwsohn/pwr


text:
using Einsum
using BenchmarkTools

srand(1234)
benchmarks = BenchmarkGroup()

### CP DECOMPOSITION BENCHMARK ###
"""
Composes a 5th-order tensor from 5 factors. This is a
CP decomposition, PARAFAC, or CANDECOMP decomposition of
tensors.
"""
function cpd(dim,r)
    A = randn(dim,r)
    B = randn(dim,r)
    C = randn(dim,r)
    D = randn(dim,r)
    E = randn(dim,r)
    @einsum X[a,b,c,d,e] := A[r,a]*B[r,b]*C[r,c]*D[r,d]*E[r,e]
    return X
end
benchmarks["cpd"] = @benchmarkable cpd(1000,10)

### COMMON FATE BENCHMARK ###
"""
Common fate bechmark proposed by @faroit
"""
function commonfate(dim,r)
    A = randn(r,dim,dim,dim)
    H = randn(r,dim)
    C = randn(r,dim)
    @einsum P[a,b,f,t,c] := A[j,a,b,f]*H[j,t]*C[j,c]
    return P
end
benchmarks["commonfate"] = @benchmarkable commonfate(30,10)


text:
using Test
import Pluto: update_save_run!, WorkspaceManager, ClientSession, ServerSession, Notebook, Cell

@testset "Workspace manager" begin
# basic functionality is already tested by the reactivity tests

    @testset "Multiple notebooks" begin

        fakeclientA = ClientSession(:fakeA, nothing)
        fakeclientB = ClientSession(:fakeB, nothing)
        🍭 = ServerSession()
        🍭.connected_clients[fakeclientA.id] = fakeclientA
        🍭.connected_clients[fakeclientB.id] = fakeclientB

        notebookA = Notebook([
            Cell("x = 3")
        ])
        fakeclientA.connected_notebook = notebookA

        notebookB = Notebook([
            Cell("x")
        ])
        fakeclientB.connected_notebook = notebookB

        @test notebookA.path != notebookB.path

        update_save_run!(🍭, notebookA, notebookA.cells[1])
        update_save_run!(🍭, notebookB, notebookB.cells[1])

        @test notebookB.cells[1].errored == true
    end

end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))

# These are the two binary objects we care about
products = Product[
    ExecutableProduct(prefix, "ffmpeg", :ffmpeg),
    ExecutableProduct(prefix, "ffprobe", :ffprobe),
    ExecutableProduct(prefix, "x264", :x264),
    ExecutableProduct(prefix, "x265", :x265),

    LibraryProduct(prefix, ["libavcodec","avcodec"], :libavcodec),
    LibraryProduct(prefix, ["libavformat","avformat"], :libavformat),
    LibraryProduct(prefix, ["libavutil","avutil"], :libavutil),
    LibraryProduct(prefix, ["libswscale","swscale"], :libswscale),
    LibraryProduct(prefix, ["libavfilter","avfilter"], :libavfilter),
    LibraryProduct(prefix, ["libavdevice","avdevice"], :libavdevice),
]

dependencies = [
    "https://github.com/JuliaPackaging/Yggdrasil/releases/download/Bzip2-v1.0.6-2/build_Bzip2.v1.0.6.jl",

    "https://github.com/ianshmean/ZlibBuilder/releases/download/v1.2.11/build_Zlib.v1.2.11.jl",
    "https://github.com/SimonDanisch/FDKBuilder/releases/download/0.1.6/build_libfdk.v0.1.6.jl",
    "https://github.com/SimonDanisch/FribidiBuilder/releases/download/0.14.0/build_fribidi.v0.14.0.jl",
    "https://github.com/JuliaGraphics/FreeTypeBuilder/releases/download/v2.9.1-4/build_FreeType2.v2.10.0.jl",
    "https://github.com/JuliaIO/LibassBuilder/releases/download/v0.14.0-2/build_libass.v0.14.0.jl",

    "https://github.com/JuliaIO/LAMEBuilder/releases/download/v3.100.0-2/build_liblame.v3.100.0.jl",

    "https://github.com/JuliaIO/OggBuilder/releases/download/v1.3.3-7/build_Ogg.v1.3.3.jl",
    "https://github.com/JuliaIO/LibVorbisBuilder/releases/download/v1.3.6-2/build_libvorbis.v1.3.6.jl",

    "https://github.com/JuliaIO/LibVPXBuilder/releases/download/v1.8.0/build_LibVPX.v1.8.0.jl",
    "https://github.com/JuliaIO/x264Builder/releases/download/v2019.5.25-static/build_x264Builder.v2019.5.25.jl",
    "https://github.com/JuliaIO/x265Builder/releases/download/v3.0.0-static/build_x265Builder.v3.0.0.jl",

    "https://github.com/JuliaIO/FFMPEGBuilder/releases/download/v4.1.0/build_FFMPEG.v4.1.0.jl"
]

for dependency in dependencies
    file = joinpath(@__DIR__, basename(dependency))
    isfile(file) || BinaryProvider.download(dependency, file)
    # it's a bit faster to run the build in an anonymous module instead of
    # starting a new julia process

    # Build the dependencies
    Mod = @eval module Anon end
    Mod.include(file)
end

# Finally, write out a deps.jl file
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
@recipe f(x::AbstractVector, gp::AbstractGP) = gp(x)
@recipe f(gp::FiniteGP) = (gp.x, gp)
@recipe function f(x::AbstractVector, gp::FiniteGP)
    length(x) == length(gp.x) ||
        throw(DimensionMismatch("length of `x` and `gp.x` has to be equal"))
    scale::Float64 = pop!(plotattributes, :ribbon_scale, 1.0)
    scale > 0.0 || error("`bandwidth` keyword argument must be non-negative")

    # compute marginals
    μ, σ2 = mean_and_var(gp)

    ribbon := scale .* sqrt.(σ2)
    fillalpha --> 0.3
    linewidth --> 2
    return x, μ
end

"""
    sampleplot([x::AbstractVector, ]f::FiniteGP; samples=1)

Plot samples from `f` versus `x` (default value: `f.x`).

Make sure to run `using Plots` before using this function.

# Example
```julia
using Plots
gp = GP(SqExponentialKernel())
sampleplot(gp(rand(5)); samples=10, markersize=5)
```
The given example plots 10 samples from the projection of the GP `gp`. The `markersize` is modified
from default of 0.5 to 5.

---
    sampleplot(x::AbstractVector, gp::AbstractGP; samples=1)

Plot samples from the finite projection `gp(x, 1e-9)` versus `x`.
"""
@userplot struct SamplePlot{X<:AbstractVector,F<:FiniteGP}
    x::X
    f::F
end

# default constructor (recipe forwards arguments as tuple)
SamplePlot((x, f)::Tuple{<:AbstractVector,<:FiniteGP}) = SamplePlot(x, f)

# `FiniteGP`s without explicit `x`
SamplePlot((f,)::Tuple{<:FiniteGP}) = SamplePlot((f.x, f))

# `AbstractGP`s with explicit `x`
# zero mean observation noise with variance 1e-9 to avoid numerical issues in
# Cholesky decomposition
SamplePlot((x, gp)::Tuple{<:AbstractVector,<:AbstractGP}) = SamplePlot((gp(x, 1e-9),))

@recipe function f(sp::SamplePlot)
    nsamples::Int = get(plotattributes, :samples, 1)
    samples = rand(sp.f, nsamples)

    # Set default attributes
    seriestype --> :line
    linealpha --> 0.2
    markershape --> :circle
    markerstrokewidth --> 0.0
    markersize --> 0.5
    markeralpha --> 0.3
    seriescolor --> "red"
    label --> ""

    return sp.x, samples
end


text:
using PGFPlots
using Base.Test

@assert success(`lualatex -v`)
using NBInclude
nbinclude(joinpath(dirname(@__FILE__), "..", "doc", "PGFPlots.ipynb"))


text:
export mask_op

function mask_op(w::judiWeights{vDT},mask::Array{mDT}) where {vDT,mDT}

    (mDT == vDT) || (mask = reshape(convert(Array{vDT},deepcopy(mask)),size(w.weights[i])))
    w_out = deepcopy(w)
    for i = 1:w.nsrc
        w_out.weights[i] = w.weights[i] .* mask
    end
    return w_out
end

function mask_op(dm::PhysicalParameter{vDT},mask::Array{mDT}) where {vDT,mDT}

    (mDT == vDT) || (mask = reshape(convert(Array{vDT},deepcopy(mask)),dm.n))
    dm_out = deepcopy(dm)
    dm_out.data = dm.data .* mask
    return dm_out

end

function mask_op(dm::Array{vDT},mask::Array{mDT}) where {vDT,mDT}

    (mDT == vDT) || (mask = convert(Array{vDT},deepcopy(mask)))
    dm_out = vec(dm) .* vec(mask)
    return dm_out

end