text:
#
# Sorting.
#

# Combinator for sorting a collection.

SortQuery(q::Query) =
    RecordQuery(q) >>
    Query(
        SortSig(false),
        Input([output(q)]),
        output(q))

SortQuery(qbase::Query, qks::Vector{Query}) =
    let q = RecordQuery(qbase >> RecordQuery(ostub(qbase), qks))
        q >>
        Query(
            SortSig(true),
            Input(domain(output(q))),
            output(qbase))
    end

SortQuery(qbase::Query, qks::Query...) =
    SortQuery(qbase, collect(qks))

# Custom ordering.

immutable SortByOrdering{O<:AbstractVector{Int}, V<:AbstractVector} <: Base.Ordering
    offs::O
    vals::V
    rev::Bool
    nullrev::Bool
end

SortByOrdering(flow::OutputFlow) =
    SortByOrdering(output(flow), column(flow))

function SortByOrdering(oty::Output, col::Column)
    @assert !isplural(oty)
    offs = offsets(col)
    vals = values(col)
    rev = decoration(oty, :rev, Bool, false)
    nullrev = decoration(oty, :nullrev, Bool, false)
    return SortByOrdering(offs, vals, rev, nullrev)
end

function Base.lt{O,V}(o::SortByOrdering{O,V}, a::Int, b::Int)
    la = o.offs[a]
    ra = o.offs[a+1]
    lb = o.offs[b]
    rb = o.offs[b+1]
    if la < ra && lb < rb
        return !o.rev ?
                isless(o.vals[la], o.vals[lb]) :
                isless(o.vals[lb], o.vals[la])
    elseif la < ra
        return o.nullrev
    elseif lb < rb
        return !o.nullrev
    else
        return false
    end
end

Base.lt{V}(o::SortByOrdering{OneTo{Int},V}, a::Int, b::Int) =
    !o.rev ?
        isless(o.vals[a], o.vals[b]) :
        isless(o.vals[b], o.vals[a])

# Sorting primitive.

immutable SortSig <: AbstractPrimitive
    haskey::Bool
end

ev(sig::SortSig, ::Input, oty::Output, iflow::InputFlow) =
    OutputFlow(oty, ev(sig, values(iflow), fields(domain(iflow))...))

function ev(sig::SortSig, dv::DataVector, oty::Output)
    if !sig.haskey
        rev = decoration(oty, :rev, Bool, false)
        return sort_impl(column(dv, 1), rev)
    else
        return sort_by_impl(column(dv, 1), fields(oty))
    end
end

function sort_impl{OPT,PLU}(col::Column{OPT,PLU}, rev::Bool)
    len = length(col)
    offs = offsets(col)
    vals = copy(values(col))
    col = Column{OPT,PLU}(offs, vals)
    cr = cursor(col)
    while !done(col, cr)
        next!(col, cr)
        sort!(cr, rev=rev)
    end
    return col
end

function sort_by_impl{OPT,PLU}(col::Column{OPT,PLU}, fs::Vector{Output})
    offs = offsets(col)
    dv = values(col)
    vals = values(dv, 1)
    perm = collect(1:length(vals))
    pcol = Column{OPT,PLU}(offs, perm)
    for k = endof(fs):-1:2
        order = SortByOrdering(fs[k], column(dv, k))
        sort_by_impl!(pcol, order)
    end
    return Column{OPT,PLU}(offs, vals[perm])
end

function sort_by_impl!(pcol::Column, order::SortByOrdering)
    cr = cursor(pcol)
    while !done(pcol, cr)
        next!(pcol, cr)
        sort!(cr, alg=MergeSort, order=order)
    end
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


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: Gerson Balbueno Bicca
#* $Id: DoublePipeIncr.mso								$
#*------------------------------------------------------------------
type DoublePipeIncr
	DoublePipeIncr()=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"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of zones",
				:Default=>2
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of incremental points",
				:Default=>3
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Options for convergence Calculations ",
				:Valid=>["Simplified","Full"],
				:Default=>"Full"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Estimated Duty",
				:Default=>70,
				:Lower=>1e-6,
				:Upper=>1e10
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight"
			)),(NComp)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Flow Direction",
				:Valid=>["counter","cocurrent"],
				:Default=>"cocurrent"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Flag for Fluid Alocation ",
				:Valid=>["outer","inner"],
				:Default=>"outer"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Inner Flow Regime ",
				:Valid=>["laminar","transition","turbulent"],
				:Default=>"laminar"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Outer Flow Regime ",
				:Valid=>["laminar","transition","turbulent"],
				:Default=>"laminar"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Laminar Flow for the Inner Side",
				:Valid=>["Hausen","Schlunder"],
				:Default=>"Hausen"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Transition Flow for the Inner Side",
				:Valid=>["Gnielinski","Hausen"],
				:Default=>"Gnielinski"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Turbulent Flow for the Inner Side",
				:Valid=>["Petukhov","SiederTate"],
				:Default=>"Petukhov"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Laminar Flow for the Outer Side",
				:Valid=>["Hausen","Schlunder"],
				:Default=>"Hausen"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Transition Flow for the OuterSide",
				:Valid=>["Gnielinski","Hausen"],
				:Default=>"Gnielinski"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Turbulent Flow for the Outer Side",
				:Valid=>["Petukhov","SiederTate"],
				:Default=>"Petukhov"
			)),
			constant(Dict{Symbol,Any}(
				:Brief=>"Pi Number",
				:Default=>3.14159265,
				:Symbol=>"\\pi"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Outside Diameter of Inner Pipe",
				:Lower=>1e-6
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Inside Diameter of Inner Pipe",
				:Lower=>1e-10
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Inside Diameter of Outer pipe",
				:Lower=>1e-10
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Effective Tube Length of one segment of Pipe",
				:Lower=>0.1,
				:Symbol=>"L_{pipe}"
			)),
			conductivity(Dict{Symbol,Any}(
				:Brief=>"Tube Wall Material Thermal Conductivity",
				:Default=>1.0,
				:Symbol=>"K_{wall}"
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Inside Fouling Resistance",
				:Unit=>"m^2*K/kW",
				:Default=>1e-6,
				:Lower=>0
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Outside Fouling Resistance",
				:Unit=>"m^2*K/kW",
				:Default=>1e-6,
				:Lower=>0
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Inner Side Outlet Nozzle Diameter",
				:Default=>0.036,
				:Lower=>10e-6
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Inner Side Inlet Nozzle Diameter",
				:Default=>0.036,
				:Lower=>10e-6
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Outer Side Outlet Nozzle Diameter",
				:Default=>0.036,
				:Lower=>10e-6
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Outer Side Inlet Nozzle Diameter",
				:Default=>0.036,
				:Lower=>10e-6
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Inner Side Inlet Nozzle Pressure Loss Coeff",
				:Default=>1.1
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Inner Side Outlet Nozzle Pressure Loss Coeff",
				:Default=>0.7
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Outer Side Inlet Nozzle Pressure Loss Coeff",
				:Default=>1.1
			)),
			positive(Dict{Symbol,Any}(
				:Brief=>"Outer Side Outlet Nozzle Pressure Loss Coeff",
				:Default=>0.7
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Inner Stream",
				:PosX=>0,
				:PosY=>0.5225,
				:Symbol=>"_{inInner}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Outer Stream",
				:PosX=>0.2805,
				:PosY=>0,
				:Symbol=>"_{inOuter}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet Inner Stream",
				:PosX=>1,
				:PosY=>0.5225,
				:Symbol=>"_{outInner}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet Outer Stream",
				:PosX=>0.7264,
				:PosY=>1,
				:Symbol=>"_{outOuter}"
			)),
			Details_Main(Dict{Symbol,Any}(
				:Brief=>"Some Details in the Heat Exchanger",
				:Symbol=>" "
			)),
			Main_DoublePipe(Dict{Symbol,Any}(
				:Brief=>"Inner Side of the Heat Exchanger",
				:Symbol=>"_{Inner}"
			)),
			Main_DoublePipe(Dict{Symbol,Any}(
				:Brief=>"Outer Side of the Heat Exchanger",
				:Symbol=>"_{Outer}"
			)),
			fill(length(Dict{Symbol,Any}(
				:Brief=>"Incremental Tube Length",
				:Symbol=>"L_{incr}"
			)),(Npoints)),
			[
				:(Outer.Properties.Average.T(1:N) = 0.5*Outer.HeatTransfer.Tlocal(1:N) + 0.5*Outer.HeatTransfer.Tlocal(2:Npoints)),
				:(Inner.Properties.Average.T(1:N) = 0.5*Inner.HeatTransfer.Tlocal(1:N) + 0.5*Inner.HeatTransfer.Tlocal(2:Npoints)),
				:(Outer.Properties.Average.P(1:N) = 0.5*Outer.PressureDrop.Plocal(1:N) + 0.5*Outer.PressureDrop.Plocal(2:Npoints)),
				:(Inner.Properties.Average.P(1:N) = 0.5*Inner.PressureDrop.Plocal(1:N) + 0.5*Inner.PressureDrop.Plocal(2:Npoints)),
				:(Inner.Properties.Wall.Twall = 0.5*Outer.Properties.Average.T + 0.5*Inner.Properties.Average.T),
				:(Outer.Properties.Wall.Twall = 0.5*Outer.Properties.Average.T + 0.5*Inner.Properties.Average.T),
				:(Outer.Properties.Average.Mw = sum(M*InletOuter.z)),
				:(Inner.Properties.Average.Mw = sum(M*InletInner.z)),
				:(Inner.Properties.Inlet.rho = PP.LiquidDensity(InletInner.T,InletInner.P,InletInner.z)),
				:(Inner.Properties.Outlet.rho = PP.LiquidDensity(OutletInner.T,OutletInner.P,OutletInner.z)),
				:(Inner.Properties.Inlet.rho = PP.VapourDensity(InletInner.T,InletInner.P,InletInner.z)),
				:(Inner.Properties.Outlet.rho = PP.VapourDensity(OutletInner.T,OutletInner.P,OutletInner.z)),
				:(Outer.Properties.Inlet.rho = PP.LiquidDensity(InletOuter.T,InletOuter.P,InletOuter.z)),
				:(Outer.Properties.Outlet.rho = PP.LiquidDensity(OutletOuter.T,OutletOuter.P,OutletOuter.z)),
				:(Outer.Properties.Inlet.rho = PP.VapourDensity(InletOuter.T,InletOuter.P,InletOuter.z)),
				:(Outer.Properties.Outlet.rho = PP.VapourDensity(OutletOuter.T,OutletOuter.P,OutletOuter.z)),
				:(Inner.Properties.Average.Cp([1:N]) = PP.LiquidCp(Inner.Properties.Average.T([1:N]),Inner.Properties.Average.P([1:N]),InletInner.z)),
				:(Inner.Properties.Average.rho([1:N]) = PP.LiquidDensity(Inner.Properties.Average.T([1:N]),Inner.Properties.Average.P([1:N]),InletInner.z)),
				:(Inner.Properties.Average.Mu([1:N]) = PP.LiquidViscosity(Inner.Properties.Average.T([1:N]),Inner.Properties.Average.P([1:N]),InletInner.z)),
				:(Inner.Properties.Average.K([1:N]) = PP.LiquidThermalConductivity(Inner.Properties.Average.T([1:N]),Inner.Properties.Average.P([1:N]),InletInner.z)),
				:(Inner.Properties.Wall.Mu([1:N]) = PP.LiquidViscosity(Inner.Properties.Wall.Twall([1:N]),Inner.Properties.Average.P([1:N]),InletInner.z)),
				:(Inner.Properties.Average.Cp([1:N]) = PP.VapourCp(Inner.Properties.Average.T([1:N]),Inner.Properties.Average.P([1:N]),InletInner.z)),
				:(Inner.Properties.Average.rho([1:N]) = PP.VapourDensity(Inner.Properties.Average.T([1:N]),Inner.Properties.Average.P([1:N]),InletInner.z)),
				:(Inner.Properties.Average.Mu([1:N]) = PP.VapourViscosity(Inner.Properties.Average.T([1:N]),Inner.Properties.Average.P([1:N]),InletInner.z)),
				:(Inner.Properties.Average.K([1:N]) = PP.VapourThermalConductivity(Inner.Properties.Average.T([1:N]),Inner.Properties.Average.P([1:N]),InletInner.z)),
				:(Inner.Properties.Wall.Mu([1:N]) = PP.VapourViscosity(Inner.Properties.Wall.Twall([1:N]),Inner.Properties.Average.P([1:N]),InletInner.z)),
				:(Outer.Properties.Average.Cp([1:N]) = PP.LiquidCp(Outer.Properties.Average.T([1:N]),Outer.Properties.Average.P([1:N]),InletOuter.z)),
				:(Outer.Properties.Average.rho([1:N]) = PP.LiquidDensity(Outer.Properties.Average.T([1:N]),Outer.Properties.Average.P([1:N]),InletOuter.z)),
				:(Outer.Properties.Average.Mu([1:N]) = PP.LiquidViscosity(Outer.Properties.Average.T([1:N]),Outer.Properties.Average.P([1:N]),InletOuter.z)),
				:(Outer.Properties.Average.K([1:N]) = PP.LiquidThermalConductivity(Outer.Properties.Average.T([1:N]),Outer.Properties.Average.P([1:N]),InletOuter.z)),
				:(Outer.Properties.Wall.Mu([1:N]) = PP.LiquidViscosity(Outer.Properties.Wall.Twall([1:N]),Outer.Properties.Average.P([1:N]),InletOuter.z)),
				:(Outer.Properties.Average.Cp([1:N]) = PP.VapourCp(Outer.Properties.Average.T([1:N]),Outer.Properties.Average.P([1:N]),InletOuter.z)),
				:(Outer.Properties.Average.rho([1:N]) = PP.VapourDensity(Outer.Properties.Average.T([1:N]),Outer.Properties.Average.P([1:N]),InletOuter.z)),
				:(Outer.Properties.Average.Mu([1:N]) = PP.VapourViscosity(Outer.Properties.Average.T([1:N]),Outer.Properties.Average.P([1:N]),InletOuter.z)),
				:(Outer.Properties.Average.K([1:N]) = PP.VapourThermalConductivity(Outer.Properties.Average.T([1:N]),Outer.Properties.Average.P([1:N]),InletOuter.z)),
				:(Outer.Properties.Wall.Mu([1:N]) = PP.VapourViscosity(Outer.Properties.Wall.Twall([1:N]),Outer.Properties.Average.P([1:N]),InletOuter.z)),
				:(Details.Q(1:N) = InletOuter.F*(Outer.HeatTransfer.Enth(1:N) - Outer.HeatTransfer.Enth(2:Npoints))),
				:(Details.Q(1:N) = InletOuter.F*(Outer.HeatTransfer.Enth(2:Npoints) - Outer.HeatTransfer.Enth(1:N))),
				:(Details.Q(1:N) = -InletInner.F*(Inner.HeatTransfer.Enth(1:N) - Inner.HeatTransfer.Enth(2:Npoints))),
				:(Details.Q(1:N) = InletInner.F*(Inner.HeatTransfer.Enth(1:N)-Inner.HeatTransfer.Enth(2:Npoints))),
				:(Details.Q(1:N) = -InletOuter.F*(Outer.HeatTransfer.Enth(1:N) - Outer.HeatTransfer.Enth(2:Npoints))),
				:(Details.Q(1:N) = -InletOuter.F*(Outer.HeatTransfer.Enth(2:Npoints) - Outer.HeatTransfer.Enth(1:N))),
				:(Inner.Properties.Inlet.Fw = sum(M*InletInner.z)*InletInner.F),
				:(Inner.Properties.Outlet.Fw = sum(M*OutletInner.z)*OutletInner.F),
				:(Outer.Properties.Inlet.Fw = sum(M*InletOuter.z)*InletOuter.F),
				:(Outer.Properties.Outlet.Fw = sum(M*OutletOuter.z)*OutletOuter.F),
				:(OutletOuter.F = InletOuter.F),
				:(OutletInner.F = InletInner.F),
				:(OutletOuter.z=InletOuter.z),
				:(OutletInner.z=InletInner.z),
				:(Details.A=Pi*DoInner*Lpipe),
				:(Lincr(1) = 0*"m"),
				:(Lincr([1:N]+1) = [1:N]*abs(Lpipe)/N),
				:(Inner.PressureDrop.fi([1:N])*Inner.PressureDrop.Re([1:N]) = 16),
				:((Inner.PressureDrop.fi([1:N])-0.0035)*(Inner.PressureDrop.Re([1:N])^0.42) = 0.264),
				:((Inner.PressureDrop.fi([1:N])-0.0035)*(Inner.PressureDrop.Re([1:N])^0.42) = 0.264),
				:(Outer.PressureDrop.fi([1:N])*Outer.PressureDrop.Re([1:N]) = 16),
				:((Outer.PressureDrop.fi([1:N])-0.0035)*(Outer.PressureDrop.Re([1:N])^0.42) = 0.264),
				:((Outer.PressureDrop.fi([1:N])-0.0035)*(Outer.PressureDrop.Re([1:N])^0.42) = 0.264),
				:(Inner.HeatTransfer.fi([1:N]) = 1/(0.79*ln(Inner.HeatTransfer.Re([1:N]))-1.64)^2),
				:(Inner.HeatTransfer.Nu([1:N]) = 3.665 + ((0.19*((DiInner/Lpipe)*Inner.HeatTransfer.Re([1:N])*Inner.HeatTransfer.PR([1:N]))^0.8)/(1+0.117*((DiInner/Lpipe)*Inner.HeatTransfer.Re([1:N])*Inner.HeatTransfer.PR([1:N]))^0.467))),
				:(Inner.HeatTransfer.Nu([1:N]) = (49.027896+4.173281*Inner.HeatTransfer.Re([1:N])*Inner.HeatTransfer.PR([1:N])*(DiInner/Lpipe))^(1/3)),
				:(Inner.HeatTransfer.fi([1:N]) = 1/(0.79*ln(Inner.HeatTransfer.Re([1:N]))-1.64)^2),
				:(Inner.HeatTransfer.Nu([1:N])*(1+(12.7*sqrt(0.125*Inner.HeatTransfer.fi([1:N]))*((Inner.HeatTransfer.PR([1:N]))^(2/3) -1))) = 0.125*Inner.HeatTransfer.fi([1:N])*(Inner.HeatTransfer.Re([1:N])-1000)*Inner.HeatTransfer.PR([1:N])),
				:(Inner.HeatTransfer.Nu([1:N]) =0.116*(Inner.HeatTransfer.Re([1:N])^(0.667)-125)*Inner.HeatTransfer.PR([1:N])^(0.333)*(1+(DiInner/Lpipe)^0.667)),
				:(Inner.HeatTransfer.fi([1:N]) = 1/(1.82*log(Inner.HeatTransfer.Re([1:N]))-1.64)^2),
				:(Inner.HeatTransfer.Nu([1:N])*(1.07+(12.7*sqrt(0.125*Inner.HeatTransfer.fi([1:N]))*((Inner.HeatTransfer.PR([1:N]))^(2/3) -1))) = 0.125*Inner.HeatTransfer.fi([1:N])*Inner.HeatTransfer.Re([1:N])*Inner.HeatTransfer.PR([1:N])),
				:(Inner.HeatTransfer.Nu([1:N]) = 0.027*(Inner.HeatTransfer.PR([1:N]))^(1/3)*(Inner.HeatTransfer.Re([1:N]))^(4/5)),
				:(Inner.HeatTransfer.fi([1:N]) = 1/(1.82*log(Inner.HeatTransfer.Re([1:N]))-1.64)^2),
				:(Outer.HeatTransfer.fi([1:N]) = 1/(0.79*ln(Outer.HeatTransfer.Re([1:N]))-1.64)^2),
				:(Outer.HeatTransfer.Nu([1:N]) = 3.665 + ((0.19*((Outer.HeatTransfer.Dh/Lpipe)*Outer.HeatTransfer.Re([1:N])*Outer.HeatTransfer.PR([1:N]))^0.8)/(1+0.117*((Outer.HeatTransfer.Dh/Lpipe)*Outer.HeatTransfer.Re([1:N])*Outer.HeatTransfer.PR([1:N]))^0.467))),
				:(Outer.HeatTransfer.Nu([1:N]) = (49.027896+4.173281*Outer.HeatTransfer.Re([1:N])*Outer.HeatTransfer.PR([1:N])*(Outer.HeatTransfer.Dh/Lpipe))^(1/3)),
				:(Outer.HeatTransfer.fi([1:N]) = 1/(0.79*ln(Outer.HeatTransfer.Re([1:N]))-1.64)^2),
				:(Outer.HeatTransfer.Nu([1:N])*(1+(12.7*sqrt(0.125*Outer.HeatTransfer.fi([1:N]))*((Outer.HeatTransfer.PR([1:N]))^(2/3) -1))) = 0.125*Outer.HeatTransfer.fi([1:N])*(Outer.HeatTransfer.Re([1:N])-1000)*Outer.HeatTransfer.PR([1:N])),
				:(Outer.HeatTransfer.Nu([1:N]) = 0.116*(Outer.HeatTransfer.Re([1:N])^(0.667)-125)*Outer.HeatTransfer.PR([1:N])^(0.333)*(1+(Outer.HeatTransfer.Dh/Lpipe)^0.667)),
				:(Outer.HeatTransfer.fi([1:N]) = 1/(0.79*ln(Outer.HeatTransfer.Re([1:N]))-1.64)^2),
				:(Outer.HeatTransfer.fi([1:N]) = 1/(1.82*log(Outer.HeatTransfer.Re([1:N]))-1.64)^2),
				:(Outer.HeatTransfer.Nu([1:N])*(1.07+(12.7*sqrt(0.125*Outer.HeatTransfer.fi([1:N]))*((Outer.HeatTransfer.PR([1:N]))^(2/3) -1))) = 0.125*Outer.HeatTransfer.fi([1:N])*Outer.HeatTransfer.Re([1:N])*Outer.HeatTransfer.PR([1:N])),
				:(Outer.HeatTransfer.Nu([1:N]) = 0.027*(Outer.HeatTransfer.PR([1:N]))^(1/3)*(Outer.HeatTransfer.Re([1:N]))^(4/5)),
				:(Outer.HeatTransfer.fi([1:N]) = 1/(1.82*log(Outer.HeatTransfer.Re([1:N]))-1.64)^2),
				:(Inner.HeatTransfer.hcoeff = (Inner.HeatTransfer.Nu*Inner.Properties.Average.K/DiInner)*Inner.HeatTransfer.Phi),
				:(Outer.HeatTransfer.hcoeff= (Outer.HeatTransfer.Nu*Outer.Properties.Average.K/Outer.HeatTransfer.Dh)*Outer.HeatTransfer.Phi),
				:(Outer.HeatTransfer.Phi = (Outer.Properties.Average.Mu/Outer.Properties.Wall.Mu)^0.14),
				:(Inner.HeatTransfer.Phi = (Inner.Properties.Average.Mu/Inner.Properties.Wall.Mu)^0.14),
				:(Outer.HeatTransfer.PR = ((Outer.Properties.Average.Cp/Outer.Properties.Average.Mw)*Outer.Properties.Average.Mu)/Outer.Properties.Average.K),
				:(Inner.HeatTransfer.PR = ((Inner.Properties.Average.Cp/Inner.Properties.Average.Mw)*Inner.Properties.Average.Mu)/Inner.Properties.Average.K),
				:(Outer.HeatTransfer.Re = (Outer.Properties.Average.rho*Outer.HeatTransfer.Vmean*Outer.HeatTransfer.Dh)/Outer.Properties.Average.Mu),
				:(Outer.PressureDrop.Re = (Outer.Properties.Average.rho*Outer.HeatTransfer.Vmean*Outer.PressureDrop.Dh)/Outer.Properties.Average.Mu),
				:(Inner.HeatTransfer.Re = (Inner.Properties.Average.rho*Inner.HeatTransfer.Vmean*Inner.HeatTransfer.Dh)/Inner.Properties.Average.Mu),
				:(Inner.PressureDrop.Re = Inner.HeatTransfer.Re),
				:(Outer.HeatTransfer.Vmean*(Outer.HeatTransfer.As*Outer.Properties.Average.rho) = Outer.Properties.Inlet.Fw),
				:(Inner.HeatTransfer.Vmean*(Inner.HeatTransfer.As*Inner.Properties.Average.rho) = Inner.Properties.Inlet.Fw),
				:(Details.Uc*((DoInner/(Inner.HeatTransfer.hcoeff*DiInner) )+(DoInner*ln(DoInner/DiInner)/(2*Kwall))+(1/(Outer.HeatTransfer.hcoeff)))=1),
				:(Details.Ud*(Rfi*(DoInner/DiInner) + Rfo + (DoInner/(Inner.HeatTransfer.hcoeff*DiInner) )+(DoInner*ln(DoInner/DiInner)/(2*Kwall))+(1/(Outer.HeatTransfer.hcoeff)))=1),
				:(Details.Qtotal = sum(Details.Q)),
				:(Details.Q = Details.Ud*Pi*DoInner*(Lpipe/N)*(Outer.Properties.Average.T - Inner.Properties.Average.T)),
				:(Details.Q = Details.Ud*Pi*DoInner*(Lpipe/N)*(Inner.Properties.Average.T - Outer.Properties.Average.T)),
				:(Details.Q = Qestimated/N),
				:(Inner.HeatTransfer.Enth([2:N]) = (1-InletInner.v)*PP.LiquidEnthalpy(Inner.HeatTransfer.Tlocal([2:N]), Inner.PressureDrop.Plocal([2:N]), InletInner.z) + InletInner.v*PP.VapourEnthalpy(Inner.HeatTransfer.Tlocal([2:N]), Inner.PressureDrop.Plocal([2:N]), InletInner.z)),
				:(Outer.HeatTransfer.Enth([2:N]) = (1-InletOuter.v)*PP.LiquidEnthalpy(Outer.HeatTransfer.Tlocal([2:N]), Outer.PressureDrop.Plocal([2:N]), InletOuter.z) + InletOuter.v*PP.VapourEnthalpy(Outer.HeatTransfer.Tlocal([2:N]), Outer.PressureDrop.Plocal([2:N]), InletOuter.z)),
				:(Inner.HeatTransfer.Enth(1) = InletInner.h),
				:(Inner.HeatTransfer.Enth(Npoints) = OutletInner.h),
				:(Inner.HeatTransfer.Tlocal(1) = InletInner.T),
				:(Inner.HeatTransfer.Tlocal(Npoints) = OutletInner.T),
				:(Inner.PressureDrop.Plocal(1) = InletInner.P),
				:(Inner.PressureDrop.Plocal(Npoints) = OutletInner.P),
				:(Outer.HeatTransfer.Enth(1) = InletOuter.h),
				:(Outer.HeatTransfer.Enth(Npoints) = OutletOuter.h),
				:(Outer.HeatTransfer.Tlocal(1) = InletOuter.T),
				:(Outer.HeatTransfer.Tlocal(Npoints) = OutletOuter.T),
				:(Outer.PressureDrop.Plocal(1) = InletOuter.P),
				:(Outer.PressureDrop.Plocal(Npoints) = OutletOuter.P),
				:(Outer.HeatTransfer.Enth(Npoints) = InletOuter.h),
				:(Outer.HeatTransfer.Enth(1) = OutletOuter.h),
				:(Outer.HeatTransfer.Tlocal(Npoints) = InletOuter.T),
				:(Outer.HeatTransfer.Tlocal(1) = OutletOuter.T),
				:(Outer.PressureDrop.Plocal(Npoints) = InletOuter.P),
				:(Outer.PressureDrop.Plocal(1) = OutletOuter.P),
				:(Outer.PressureDrop.Pdrop = Outer.PressureDrop.Pd_fric(Npoints)+Outer.PressureDrop.Pdnozzle_in+Outer.PressureDrop.Pdnozzle_out),
				:(Outer.PressureDrop.Pd_fric(2:Npoints) = (2*Outer.PressureDrop.fi*Lincr(2:Npoints)*Outer.Properties.Average.rho*Outer.HeatTransfer.Vmean^2)/(Outer.PressureDrop.Dh*Outer.HeatTransfer.Phi)),
				:(Outer.PressureDrop.Pd_fric(1) = 0*"kPa"),
				:(Outer.PressureDrop.Plocal([1:N]+1) = Outer.PressureDrop.Plocal(1) - Outer.PressureDrop.Pd_fric([1:N]+1)),
				:(Outer.PressureDrop.Pdrop = Outer.PressureDrop.Pd_fric(1)+Outer.PressureDrop.Pdnozzle_in+Outer.PressureDrop.Pdnozzle_out),
				:(Outer.PressureDrop.Pd_fric([1:N]) = (2*Outer.PressureDrop.fi([1:N])*Lincr(1+N-[1:N])*Outer.Properties.Average.rho([1:N])*Outer.HeatTransfer.Vmean([1:N])^2)/(Outer.PressureDrop.Dh*Outer.HeatTransfer.Phi([1:N]))),
				:(Outer.PressureDrop.Pd_fric(Npoints) = 0*"kPa"),
				:(Outer.PressureDrop.Plocal([1:N]) = Outer.PressureDrop.Plocal(Npoints) - Outer.PressureDrop.Pd_fric([1:N]+1)),
				:(Inner.PressureDrop.Pdrop = Inner.PressureDrop.Pd_fric(Npoints)+Inner.PressureDrop.Pdnozzle_in+Inner.PressureDrop.Pdnozzle_out),
				:(Inner.PressureDrop.Pd_fric(2:Npoints) = (2*Inner.PressureDrop.fi*Lincr(2:Npoints)*Inner.Properties.Average.rho*Inner.HeatTransfer.Vmean^2)/(DiInner*Inner.HeatTransfer.Phi)),
				:(Inner.PressureDrop.Pd_fric(1) = 0*"kPa"),
				:(Inner.PressureDrop.Plocal([1:N]+1) = Inner.PressureDrop.Plocal(1) - Inner.PressureDrop.Pd_fric([1:N]+1)),
				:(Inner.PressureDrop.Vnozzle_in = Inner.Properties.Inlet.Fw/(Inner.Properties.Inlet.rho*(0.25*Pi*Dinozzle_Inner^2))),
				:(Inner.PressureDrop.Vnozzle_out = Inner.Properties.Outlet.Fw/(Inner.Properties.Outlet.rho*(0.25*Pi*Donozzle_Inner^2))),
				:(Outer.PressureDrop.Vnozzle_in = Outer.Properties.Inlet.Fw/(Outer.Properties.Inlet.rho*(0.25*Pi*Dinozzle_Outer^2))),
				:(Outer.PressureDrop.Vnozzle_out = Outer.Properties.Outlet.Fw/(Outer.Properties.Outlet.rho*(0.25*Pi*Donozzle_Outer^2))),
				:(Inner.PressureDrop.Pdnozzle_in = 0.5*InnerKinlet*Inner.Properties.Inlet.rho*Inner.PressureDrop.Vnozzle_in^2),
				:(Inner.PressureDrop.Pdnozzle_out = 0.5*InnerKoutlet*Inner.Properties.Outlet.rho*Inner.PressureDrop.Vnozzle_out^2),
				:(Outer.PressureDrop.Pdnozzle_in = 0.5*OuterKinlet*Outer.Properties.Inlet.rho*Outer.PressureDrop.Vnozzle_in^2),
				:(Outer.PressureDrop.Pdnozzle_out = 0.5*OuterKoutlet*Outer.Properties.Outlet.rho*Outer.PressureDrop.Vnozzle_out^2),
				:(Inner.PressureDrop.RVsquare_in = Inner.Properties.Inlet.rho*(Inner.PressureDrop.Vnozzle_in)^2),
				:(Inner.PressureDrop.RVsquare_out = Inner.Properties.Outlet.rho*(Inner.PressureDrop.Vnozzle_out)^2),
				:(Outer.PressureDrop.RVsquare_in = Outer.Properties.Inlet.rho*(Outer.PressureDrop.Vnozzle_in)^2),
				:(Outer.PressureDrop.RVsquare_out = Outer.Properties.Outlet.rho*(Outer.PressureDrop.Vnozzle_out)^2),
			],
			[
				"Outer	Stream Average Temperature","Inner Stream Average Temperature","Outer Stream Average Pressure","Inner Stream Average Pressure","Inner Stream Wall Temperature","Outer Stream Wall Temperature","Outer Stream Average Molecular Weight","Inner Stream Average Molecular Weight","Inlet Mass Density Inner Stream","Outlet Mass Density Inner Stream","Inlet Mass Density Inner Stream","Outlet Mass Density Inner Stream","Inlet Mass Density Outer Stream","Outlet Mass Density Outer Stream","Inlet Mass Density Outer Stream","Outlet Mass Density Outer Stream","Average Heat Capacity Inner Stream","Average Mass Density Inner Stream","Average Viscosity Inner Stream","Average	Conductivity Inner Stream","Viscosity Inner Stream at wall temperature","Average Heat Capacity InnerStream","Average Mass Density Inner Stream","Average Viscosity Inner Stream","Average Conductivity Inner Stream","Viscosity Inner Stream at wall temperature","Average Heat Capacity Outer Stream","Average Mass Density Outer Stream","Average Viscosity Outer Stream","Average Conductivity Outer Stream","Viscosity Outer Stream at wall temperature","Average Heat Capacity Outer Stream","Average Mass Density Outer Stream","Average Viscosity Outer Stream","Average Conductivity Outer Stream","Viscosity Outer Stream at wall temperature","Energy Balance Outer Stream in cocurrent flow","Energy Balance Outer Stream in counter flow","Energy Balance Inner Stream","Energy Balance Hot Stream","Energy Balance Cold Stream in cocurrent flow","Energy Balance Cold Stream in counter flow","Flow Mass Inlet Inner Stream","Flow Mass Outlet Inner Stream","Flow Mass Inlet Outer Stream","Flow Mass Outlet Outer Stream","Molar Balance Outer Stream","Molar Balance Inner Stream","Outer Stream Molar Fraction Constraint","InnerStream Molar Fraction Constraint","Total Exchange Surface Area for one segment of pipe","Pipe Initial Length from Left to Right - OBS: Left: Always Inlet inner side","Incremental Length","Inner Side Friction Factor for Pressure Drop - laminar Flow","using Turbulent Flow - to be implemented","Inner Side Friction Factor - Turbulent Flow","Outer Side Friction Factor - laminar Flow","using Turbulent Flow - Transition Flow must be implemented","Outer Side Friction Factor - Turbulent Flow","Inner Side Friction Factor for Heat Transfer - laminar Flow","Nusselt Number","Nusselt Number","Inner Side Friction Factor for Heat Transfer - transition Flow","Nusselt Number","Nusselt Number","Inner Side Friction Factor for Heat Transfer - turbulent Flow","Nusselt Number","Nusselt Number","Inner Side Friction Factor for Heat Transfer - turbulent Flow","Outer Side Friction Factor for Heat Transfer - laminar Flow","Nusselt Number","Nusselt Number","Outer Side Friction Factor for Heat Transfer - transition Flow","Nusselt Number","Nusselt Number","Outer Side Friction Factor for Heat Transfer - transition Flow","Outer Side Friction Factor for Heat Transfer - turbulent Flow","Nusselt Number","Nusselt Number","Outer Side Friction Factor for Heat Transfer - turbulent Flow","Inner Pipe Film Coefficient","Outer Pipe Film Coefficient","Outer Pipe Phi correction","Inner Pipe Phi correction","Outer Pipe Prandtl Number","Inner Pipe Prandtl Number","Outer Pipe Reynolds Number for Heat Transfer","Outer Pipe Reynolds Number for Pressure Drop","Inner Pipe Reynolds Number for Heat Transfer","Inner Pipe Reynolds Number for Pressure Drop","Outer Pipe Velocity","Inner Pipe Velocity","Overall Heat Transfer Coefficient Clean","Overall Heat Transfer Coefficient Dirty","Total Duty","Incremental Duty","Incremental Duty","Duty","Incremental Enthalpy Inner Stream","Incremental Enthalpy Outer Stream","Enthalpy of Inner Side - Inlet Boundary","Enthalpy of inner Side - Outlet Boundary","Temperature of Inner Side - Inlet Boundary","Temperature of Inner Side - Outlet Boundary","Pressure of Inner Side - Inlet Boundary","Pressure of Inner Side - Outlet Boundary","Enthalpy of Outer Side - Inlet Boundary","Enthalpy of Outer Side - Outlet Boundary","Temperature of Outer Side - Inlet Boundary","Temperature of Outer Side - Outlet Boundary","Pressure of Outer Side - Inlet Boundary","Pressure of Outer Side - Outlet Boundary","Enthalpy of Outer Side - Inlet Boundary","Enthalpy of Outer Side - Outlet Boundary","Temperature of Outer Side - Inlet Boundary","Temperature of Outer Side - Outlet Boundary","Pressure of Outer Side - Inlet Boundary","Pressure of Outer Side - Outlet Boundary","Total Pressure Drop Outer Stream","Outer Pipe Pressure Drop for friction","Outer Pipe Pressure Drop for friction","Outer Pipe Local Pressure","Total Pressure Drop Outer Stream","Outer Pipe Pressure Drop for friction","Outer Pipe Pressure Drop for friction","Outer Pipe Local Pressure","Total Pressure Drop Inner Stream","Inner Pipe Pressure Drop for friction","Inner Pipe Pressure Drop for friction","Inner Pipe Local Pressure","Velocity Inner Side Inlet Nozzle","Velocity Inner Side Outlet Nozzle","Velocity Outer Side Inlet Nozzle","Velocity Outer Side Outlet Nozzle","Pressure Drop Inner Side Inlet Nozzle","Pressure Drop Inner Side Outlet Nozzle","Pressure Drop Outer Side Inlet Nozzle","Pressure Drop Outer Side Outlet Nozzle","Inner Side Inlet Nozzle rho-V^2","Inner Side Outlet Nozzle rho-V^2","Outer Side Inlet Nozzle rho-V^2","Outer Side Outlet Nozzle rho-V^2",
			],
			[:PP,:NComp,:N,:Npoints,:CalculationApproach,:Qestimated,:M,:FlowDirection,:HotSide,:innerFlowRegime,:outerFlowRegime,:InnerLaminarCorrelation,:InnerTransitionCorrelation,:InnerTurbulentCorrelation,:OuterLaminarCorrelation,:OuterTransitionCorrelation,:OuterTurbulentCorrelation,:Pi,:DoInner,:DiInner,:DiOuter,:Lpipe,:Kwall,:Rfi,:Rfo,:Donozzle_Inner,:Dinozzle_Inner,:Donozzle_Outer,:Dinozzle_Outer,:InnerKinlet,:InnerKoutlet,:OuterKinlet,:OuterKoutlet,],
			[:InletInner,:InletOuter,:OutletInner,:OutletOuter,:Details,:Inner,:Outer,:Lincr,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	N::DanaInteger
	Npoints::DanaInteger
	CalculationApproach::DanaSwitcher
	Qestimated::power
	M::Array{molweight}
	FlowDirection::DanaSwitcher
	HotSide::DanaSwitcher
	innerFlowRegime::DanaSwitcher
	outerFlowRegime::DanaSwitcher
	InnerLaminarCorrelation::DanaSwitcher
	InnerTransitionCorrelation::DanaSwitcher
	InnerTurbulentCorrelation::DanaSwitcher
	OuterLaminarCorrelation::DanaSwitcher
	OuterTransitionCorrelation::DanaSwitcher
	OuterTurbulentCorrelation::DanaSwitcher
	Pi::constant
	DoInner::length
	DiInner::length
	DiOuter::length
	Lpipe::length
	Kwall::conductivity
	Rfi::positive
	Rfo::positive
	Donozzle_Inner::length
	Dinozzle_Inner::length
	Donozzle_Outer::length
	Dinozzle_Outer::length
	InnerKinlet::positive
	InnerKoutlet::positive
	OuterKinlet::positive
	OuterKoutlet::positive
	InletInner::stream
	InletOuter::stream
	OutletInner::streamPH
	OutletOuter::streamPH
	Details::Details_Main
	Inner::Main_DoublePipe
	Outer::Main_DoublePipe
	Lincr::Array{length}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export DoublePipeIncr
function set(in::DoublePipeIncr)
	#"Number of incremental points"
	Npoints = N+1
	 #"Component Molecular Weight"
	M = PP.MolecularWeight()
	 #"Pi Number"
	Pi = 3.14159265
	 #"Inner Pipe Cross Sectional Area for Flow"
	Inner.HeatTransfer.As=0.25*Pi*DiInner*DiInner
	 #"Outer Pipe Cross Sectional Area for Flow"
	Outer.HeatTransfer.As=0.25*Pi*(DiOuter*DiOuter - DoInner*DoInner)
	 #"Inner Pipe Hydraulic Diameter for Heat Transfer"
	Inner.HeatTransfer.Dh=DiInner
	 #"Outer Pipe Hydraulic Diameter for Heat Transfer"
	Outer.HeatTransfer.Dh=(DiOuter*DiOuter-DoInner*DoInner)/DoInner
	 #"Inner Pipe Hydraulic Diameter for Pressure Drop"
	Inner.PressureDrop.Dh=DiInner
	 #"Outer Pipe Hydraulic Diameter for Pressure Drop"
	Outer.PressureDrop.Dh=DiOuter-DoInner

end
function setEquationFlow(in::DoublePipeIncr)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	if InletInner.v == 0
		addEquation(9)
		addEquation(10)
	else
		addEquation(11)
		addEquation(12)
	end
	if InletOuter.v == 0
		addEquation(13)
		addEquation(14)
	else
		addEquation(15)
		addEquation(16)
	end
	if InletInner.v == 0
		addEquation(17)
		addEquation(18)
		addEquation(19)
		addEquation(20)
		addEquation(21)
	else
		addEquation(22)
		addEquation(23)
		addEquation(24)
		addEquation(25)
		addEquation(26)
	end
	if InletOuter.v == 0
		addEquation(27)
		addEquation(28)
		addEquation(29)
		addEquation(30)
		addEquation(31)
	else
		addEquation(32)
		addEquation(33)
		addEquation(34)
		addEquation(35)
		addEquation(36)
	end
	let switch=HotSide
		if HotSide==InletInner.T > InletOuter.T
			set(switch,"inner")
		end
		if HotSide==InletInner.T < InletOuter.T
			set(switch,"outer")
		end
		if switch=="outer"
			let switch=FlowDirection
				if switch=="cocurrent"
					addEquation(37)
				elseif switch=="counter"
					addEquation(38)
				end
			end
			addEquation(39)
		elseif switch=="inner"
			addEquation(40)
			let switch=FlowDirection
				if switch=="cocurrent"
					addEquation(41)
				elseif switch=="counter"
					addEquation(42)
				end
			end
		end
	end
	addEquation(43)
	addEquation(44)
	addEquation(45)
	addEquation(46)
	addEquation(47)
	addEquation(48)
	addEquation(49)
	addEquation(50)
	addEquation(51)
	addEquation(52)
	addEquation(53)
	let switch=innerFlowRegime
		if innerFlowRegime==Inner.PressureDrop.Re(i) > 2300
			set(switch,"transition")
		end
		if innerFlowRegime==Inner.PressureDrop.Re(i) < 2300
			set(switch,"laminar")
		end
		if innerFlowRegime==Inner.PressureDrop.Re(i) > 10000
			set(switch,"turbulent")
		end
		if innerFlowRegime==Inner.PressureDrop.Re(i) < 10000
			set(switch,"transition")
		end
		if switch=="laminar"
			addEquation(54)
		elseif switch=="transition"
			addEquation(55)
		elseif switch=="turbulent"
			addEquation(56)
		end
	end
	let switch=outerFlowRegime
		if outerFlowRegime==Outer.PressureDrop.Re(i) > 2300
			set(switch,"transition")
		end
		if outerFlowRegime==Outer.PressureDrop.Re(i) < 2300
			set(switch,"laminar")
		end
		if outerFlowRegime==Outer.PressureDrop.Re(i) > 10000
			set(switch,"turbulent")
		end
		if outerFlowRegime==Outer.PressureDrop.Re(i) < 10000
			set(switch,"transition")
		end
		if switch=="laminar"
			addEquation(57)
		elseif switch=="transition"
			addEquation(58)
		elseif switch=="turbulent"
			addEquation(59)
		end
	end
	let switch=innerFlowRegime
		if innerFlowRegime==Inner.HeatTransfer.Re(i) > 2300
			set(switch,"transition")
		end
		if innerFlowRegime==Inner.HeatTransfer.Re(i) < 2300
			set(switch,"laminar")
		end
		if innerFlowRegime==Inner.HeatTransfer.Re(i) > 10000
			set(switch,"turbulent")
		end
		if innerFlowRegime==Inner.HeatTransfer.Re(i) < 10000
			set(switch,"transition")
		end
		if switch=="laminar"
			addEquation(60)
			let switch=InnerLaminarCorrelation
				if switch=="Hausen"
					addEquation(61)
				elseif switch=="Schlunder"
					addEquation(62)
				end
			end
		elseif switch=="transition"
			addEquation(63)
			let switch=InnerTransitionCorrelation
				if switch=="Gnielinski"
					addEquation(64)
				elseif switch=="Hausen"
					addEquation(65)
				end
			end
		elseif switch=="turbulent"
			let switch=InnerTurbulentCorrelation
				if switch=="Petukhov"
					addEquation(66)
					addEquation(67)
				elseif switch=="SiederTate"
					addEquation(68)
					addEquation(69)
				end
			end
		end
	end
	let switch=outerFlowRegime
		if outerFlowRegime==Outer.HeatTransfer.Re(i) > 2300
			set(switch,"transition")
		end
		if outerFlowRegime==Outer.HeatTransfer.Re(i) < 2300
			set(switch,"laminar")
		end
		if outerFlowRegime==Outer.HeatTransfer.Re(i) > 10000
			set(switch,"turbulent")
		end
		if outerFlowRegime==Outer.HeatTransfer.Re(i) < 10000
			set(switch,"transition")
		end
		if switch=="laminar"
			addEquation(70)
			let switch=OuterLaminarCorrelation
				if switch=="Hausen"
					addEquation(71)
				elseif switch=="Schlunder"
					addEquation(72)
				end
			end
		elseif switch=="transition"
			let switch=OuterTransitionCorrelation
				if switch=="Gnielinski"
					addEquation(73)
					addEquation(74)
				elseif switch=="Hausen"
					addEquation(75)
					addEquation(76)
				end
			end
		elseif switch=="turbulent"
			let switch=OuterTurbulentCorrelation
				if switch=="Petukhov"
					addEquation(77)
					addEquation(78)
				elseif switch=="SiederTate"
					addEquation(79)
					addEquation(80)
				end
			end
		end
	end
	addEquation(81)
	addEquation(82)
	addEquation(83)
	addEquation(84)
	addEquation(85)
	addEquation(86)
	addEquation(87)
	addEquation(88)
	addEquation(89)
	addEquation(90)
	addEquation(91)
	addEquation(92)
	addEquation(93)
	addEquation(94)
	addEquation(95)
	let switch=CalculationApproach
		if switch=="Full"
			let switch=HotSide
				if HotSide==InletInner.T > InletOuter.T
					set(switch,"inner")
				end
				if HotSide==InletInner.T < InletOuter.T
					set(switch,"outer")
				end
				if switch=="outer"
					addEquation(96)
				elseif switch=="inner"
					addEquation(97)
				end
			end
		elseif switch=="Simplified"
			addEquation(98)
		end
	end
	addEquation(99)
	addEquation(100)
	addEquation(101)
	addEquation(102)
	addEquation(103)
	addEquation(104)
	addEquation(105)
	addEquation(106)
	let switch=FlowDirection
		if switch=="cocurrent"
			addEquation(107)
			addEquation(108)
			addEquation(109)
			addEquation(110)
			addEquation(111)
			addEquation(112)
		elseif switch=="counter"
			addEquation(113)
			addEquation(114)
			addEquation(115)
			addEquation(116)
			addEquation(117)
			addEquation(118)
		end
	end
	let switch=FlowDirection
		if switch=="cocurrent"
			addEquation(119)
			addEquation(120)
			addEquation(121)
			# FIXME: NOZZLE PRESSURE DROP MUST BE ADDED
			addEquation(122)
		elseif switch=="counter"
			addEquation(123)
			addEquation(124)
			addEquation(125)
			# FIXME: NOZZLE PRESSURE DROP MUST BE ADDED
			addEquation(126)
		end
	end
	addEquation(127)
	addEquation(128)
	addEquation(129)
	# FIXME: NOZZLE PRESSURE DROP MUST BE ADDED
	addEquation(130)
	addEquation(131)
	addEquation(132)
	addEquation(133)
	addEquation(134)
	addEquation(135)
	addEquation(136)
	addEquation(137)
	addEquation(138)
	addEquation(139)
	addEquation(140)
	addEquation(141)
	addEquation(142)
end
function atributes(in::DoublePipeIncr,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/DoublePipe"
	fields[:Brief]="Incremental Double Pipe Heat Exchanger. "
	fields[:Info]="Incremental approach for a single double pipe heat exchanger. "
	drive!(fields,_)
	return fields
end
DoublePipeIncr(_::Dict{Symbol,Any})=begin
	newModel=DoublePipeIncr()
	newModel.attributes=atributes(newModel,_)
	newModel
end

ext: jl


text:
@non_differentiable id(V::VectorSpace)
@non_differentiable isomorphism(cod::VectorSpace, dom::VectorSpace)
@non_differentiable domain(W::TensorKit.HomSpace)
@non_differentiable domain(t::AbstractTensorMap)
@non_differentiable domain(t::AbstractTensorMap, i)
@non_differentiable codomain(W::TensorKit.HomSpace)
@non_differentiable codomain(t::AbstractTensorMap)
@non_differentiable codomain(t::AbstractTensorMap, i)

convert_to_array(X::AbstractTensorMap) = convert(Array, X)

function rrule(::typeof(convert_to_array), X::AbstractTensorMap)

    fwd = convert_to_array(X)
    space_X = space(X)

    function convert_to_array_pushback(f̄wd)
        return NoTangent(), TensorMap(f̄wd, space_X)
    end

    return fwd, convert_to_array_pushback
end

"""
    rrule(::typeof(TensorMap), data::A, codom::ProductSpace, dom::ProductSpace) where A <: DenseMatrix

    rrule for TensorMap construction. only did the simplest case.
    FIXME. seems buggy
"""
function rrule(::typeof(TensorMap), data::A, codom::ProductSpace, dom::ProductSpace) where A <: DenseMatrix
    fwd = TensorMap(data, codom, dom)
    function TensorMap_pushback(f̄wd)
        return NoTangent(), f̄wd.data, NoTangent(), NoTangent()
    end
    return fwd, TensorMap_pushback
end

"""
    rrule(::typeof(tr), A::TensorMap)

    rrule for trace operation.
"""
function rrule(::typeof(tr), A::TensorMap)
    fwd = tr(A)
    function tr_pushback(f̄wd)
        Ā = f̄wd * id(domain(A))
        return NoTangent(), Ā
    end
    return fwd, tr_pushback
end

"""
    elem_mult(a::AbstractTensorMap,b::AbstractTensorMap)

    Element-wise multiplication between TensorMaps. Copied from TensorKitAD.jl.
"""
function elem_mult(a::AbstractTensorMap,b::AbstractTensorMap)
    dst = similar(a);
    for (k,block) in blocks(dst)
        copyto!(block,blocks(a)[k].*blocks(b)[k]);
    end
    dst
end

"""
    rrule(::typeof(elem_mult), a::AbstractTensorMap, b::AbstractTensorMap)
"""
function rrule(::typeof(elem_mult), a::AbstractTensorMap, b::AbstractTensorMap)
    fwd = similar(a);
    for (k,block) in blocks(fwd)
        copyto!(block,blocks(a)[k].*blocks(b)[k]);
    end

    function elem_mult_pushback(f̄wd)
        ā, b̄ = similar(a), similar(b)
        for (k, block) in blocks(f̄wd)
            copyto!(blocks(ā)[k], block .* conj.(blocks(b)[k]))
            copyto!(blocks(b̄)[k], block .* conj.(blocks(a)[k]))
        end

        return NoTangent(), ā, b̄
    end
    return fwd, elem_mult_pushback
end


text:
using GLM
model = lm(@formula(r ~ Protocol + Wall),comp)
nullmodel = lm(@formula(r ~ Protocol ),comp)
F1 = ftest(nullmodel.model, model.model)

model = lm(@formula(r ~ Protocol + Wall),comp)
nullmodel = lm(@formula(r ~ Wall ),comp)
ftest(nullmodel.model, model.model)

model = lm(@formula(T ~ Protocol + Wall),comp)
nullmodel = lm(@formula(T ~ Protocol ),comp)
ftest(nullmodel.model, model.model)

model = lm(@formula(T ~ Protocol + Wall),comp)
nullmodel = lm(@formula(T ~ Wall ),comp)
ftest(nullmodel.model, model.model)

###

model = lm(@formula(lapse ~ Protocol + Wall),comp)
nullmodel = lm(@formula(lapse ~ Protocol ),comp)
ftest(nullmodel.model, model.model)

model = lm(@formula(lapse ~ Protocol + Wall),comp)
nullmodel = lm(@formula(lapse ~ Wall ),comp)
ftest(nullmodel.model, model.model)
##
Rate_model = lm(@formula(r ~ 1 + Protocol + Wall),comp);
Rate_nullmodel = lm(@formula(r ~ 1 ),comp);
Rate_F = ftest(Rate_nullmodel.model, Rate_model.model);
##
Tresh_model = lm(@formula(T ~ 1 + Protocol + Wall),comp);
Tresh_nullmodel = lm(@formula(T ~ 1 ),comp);
Tresh_F = ftest(Tresh_nullmodel.model, Tresh_model.model);
##
Lapse_model = lm(@formula(lapse ~ 1 + Protocol + Wall),comp);
Lapse_nullmodel = lm(@formula(lapse ~ 1 ),comp);
Lapse_F = ftest(Lapse_nullmodel.model, Lapse_model.model);
##
args, kwargs = series2D(
    Recombinase.discrete(Recombinase.prediction),
    comp,
    Recombinase.Group(color = :Wall),
    select = (:Protocol, :r)
    )
scatter(args...;kwargs...,
    annotate =(0.8,3,text("F-test: p < 10^-7",6,:black,:center)),
    grid = false,
    linewidth = 2.5,
    markersize = 8,
    thickness_scaling = 2.5,
    tick_orientation = :out,
    size = (809,500),
    legend = false)
#scatter!(args...;kwargs...,markersize = 8,legend = false)
##
savefig("/Volumes/GoogleDrive/My Drive/Reports for Zach/NewProtocols/Fitting/rate_modulation.pdf")
##
figures = "/Users/dariosarra/Documents/Lab/Mainen/Presentations/Walton_visit/figures/unprocessed_plots"
savefig(joinpath(figures,"Rate_scatter.pdf"))
##
args, kwargs = series2D(
    Recombinase.discrete(Recombinase.prediction),
    comp,
    Recombinase.Group(color = :Wall),
    select = (:Protocol, :T)
    )
scatter(args...;kwargs...,
    annotate =(0.8,8.7,text("F-test: p = 0.3991",6,:black,:center)),
    grid = false,
    linewidth = 2.5,
    markersize = 8,
    thickness_scaling = 2.5,
    tick_orientation = :out,
    size = (809,500),
    legend = false)
# scatter!(args...;kwargs...,markersize = 8,legend = false)
##
figures = "/Users/dariosarra/Documents/Lab/Mainen/Presentations/Walton_visit/figures/unprocessed_plots"
savefig(joinpath(figures,"Threshold_scatter.pdf"))
##
savefig("/Volumes/GoogleDrive/My Drive/Reports for Zach/NewProtocols/Fitting/tresh_modulation.pdf")
##
args, kwargs = series2D(
    Recombinase.discrete(Recombinase.prediction),
    comp,
    Recombinase.Group(color = :Wall),
    select = (:Protocol, :lapse)
    )
scatter(args...;kwargs...,
    annotate =(0.8,0.095,text("F-test: p = 0.9067",6,:black,:center)),
    grid = false,
    linewidth = 2.5,
    markersize = 8,
    thickness_scaling = 2.5,
    tick_orientation = :out,
    size = (809,500),
    legend = false)
#scatter!(args...;kwargs...,markersize = 8,legend = false)
##
figures = "/Users/dariosarra/Documents/Lab/Mainen/Presentations/Walton_visit/figures/unprocessed_plots"
savefig(joinpath(figures,"Lapse_scatter.pdf"))
##
savefig("/Volumes/GoogleDrive/My Drive/Reports for Zach/NewProtocols/Fitting/lapse_modulation.pdf")


text:
using Mimi
using MimiRICE2010

m = get_model()
run(m)

explore(m)


text:
"This script is intended to manipulate the multidimensional arrays for the image transformation"

import Images

export imgload, chsplit

"`imgload(path_string)` returns the array of values of a monochromatic or polychromatic image. In the first case is onlye an N×N matrix, in the second case is a collection of three N×N matrices."
function imgload(path::String)
    return Images.load(path)
end

"`chsplit(img, c)` returns the c-channel of and image _img_. It requires a polychromatic imported image from the `imgload` function, and a channel c ∈ {1, 2, 3} whose equivalences are 1 = r (red), 2 = g (green), 3 = b (blue)."
function chsplit(img, c)
    chs = Images.channelview(img)
    if c == 1
        return chs[1,:,:]
    elseif c == 2
        return chs[2,:,:]
    elseif c == 3
        return chs[3,:,:]
    end
end

"`imgrep(r, g, b)` takes the channel arrays (original or transformated), and returns a stacked representation of dimensions N×N×3, suitable for plotting in the `imshow` command of PyPlot"
function imgrep(r,g,b)
    return permutedims(Images.StackedView(r, g, b), [3, 2, 1])
end


text:
len(n) = length(string(n))

c = 0
n = 1
while true
  global c, n
  if len(BigInt(9^n)) < n
    break
  end
  for i = 0:9
    maybe = BigInt(i^n)
    if len(maybe) == n
      println(i, "^", n, "=", maybe)
      c += 1
    end
  end
  n += 1
end
println(c)

max_stars_repo_name: pdbartlett/misc-stuff

max_issues_repo_name: pdbartlett/misc-stuff

max_forks_repo_name: pdbartlett/misc-stuff


text:
function head(tree :: TreeNode)
    if isterminal(tree)
        tree
    else
        nodecat = features(tree.data[1].cat)
        headfeat = catfeature(nodecat)
        #head(filter(t -> catfeature(t.data[1].cat[1]) == headfeat, tree.children)[1])
        if selects_right(features(tree.children[1].data[1].cat)[1])
            head(tree.children[1])
        else
            head(tree.children[2])
        end
    end
end

"""
This function assumes a binary MG derivation tree.
Give a list of dependents of the node at the root of tree. Dependents are represented
as a list of indices into the leaves (indices of their linear order).
"""
function dependents(tree :: TreeNode)
    nodecat = features(tree.data[1].cat)
    if length(nodecat) > 1 # has selector features
        dep = if selects_right(nodecat[1])
            head(tree.parent.children[2])
        else
            head(tree.parent.children[1])
        end
        [dependents(tree.parent); dep]
    else
        TreeNode[]
    end
end

"""
Given a tree, gives an unlabelled dependency graph over the words in the tree.
Return type is SimpleDiGraph.
"""
function dependencies{T}(tree :: TreeNode{T};
                         closed_under_transitivity=false,
                         directed=true)
    non_null_head(t) = !isempty(t.data[2][2])
    leaves = filter(non_null_head, leafs(tree))
    g = if directed
        DiGraph(length(leaves))
    else
        Graph(length(leaves))
    end

    for (i, leaf) in enumerate(leaves)
        dents = dependents(leaf)
        for (j, dent) in enumerate(leaves)
            if dent in dents
                add_edge!(g, i, j)
            end
        end
    end

    if closed_under_transitivity
        if directed
            transitiveclosure(g)
        else
            # transitiveclosure only defined for di-graphs
            Graph(transitiveclosure(DiGraph(g)))
        end
    else g end
end

"""
If a dependency matrix has dependencies (i,j) where j-i > 2, and doesn't have dependencies inside i,j that give a binary tree, add them.

It works like this:
Let (i,j) be a cell in the matrix. If i - j > 3 then (i,j) spans more than 2 words. We check if either cells (i+1, j) or (i, j-1) are true. If not, then (i,j) is not binarized, and set both of those cells to true.
"""
function binarize_dependency_matrix(dm :: AbstractMatrix{Bool})
    n = size(dm,1)
    new_dm = copy(dm)
    # helper methods
    # is_binarized checks if i,j is binarized, uses bm for memoization
    is_binarized(i, j) = any([dm[i, j-1], dm[i+1, j]])
    # lt compares ranges (i,j) by their span length j-i
    lt(x1, x2) = x1[2]-x1[1] < x2[2]-x2[1]
    # push_sorted pushes an item into a sorted list using our lt
    push_sorted!(l, x) = insert!(l, searchsortedfirst(l, x, lt=lt), x)
    # get indices of dependencies spanning more than 2 words & sort it
    long_dependencies = [(i,j) for i in 1:n for j in i+3:n if dm[i,j]]
    sort!(long_dependencies, lt=lt)
    # binarize
    while !isempty(long_dependencies)
        i, j = pop!(long_dependencies)
        if !is_binarized(i, j)
            # add dependencies
            new_dm[i, j-1] = true
            new_dm[i+1, j] = true
            if j-i > 3
                push_sorted!(long_dependencies, (i, j-1))
                push_sorted!(long_dependencies, (i+1, j))
            end
        end
    end
    new_dm
end

"Converts a dependency graph with n vertices to an (n+1)×(n+1) dependency matrix"
function graph_to_dependency_matrix(g :: DiGraph{Int64})
    n = nv(g) + 1
    m = zeros(Matrix{Bool}(n, n))
    function dependencies_span(v)
        out_spans = map(dependencies_span, outneighbors(g, v))
        is = [s[1] for s in out_spans]
        js = [s[2] for s in out_spans]
        i = minimum([v  ; outneighbors(g, v)               ; is])
        j = maximum([v+1; [w+1 for w in outneighbors(g, v)]; js])
        i, j
    end
    for v in vertices(g)
        i, j = dependencies_span(v)
        m[i, j]   = true
        m[v, v+1] = true
    end
    m
end


text:
#
# Returns a view of a coordinate vector corresponding to the atoms
# of a molecule with index i. n is the number of atoms of the molecule
#
function viewmol(i :: Int, x :: Vector{T}, n :: Int) where T
  first = (i-1)* n + 1
  last = first + n - 1
  return @view(x[first:last])
end

# From the structure of the solute or solvent
viewmol(i :: Int, x :: Vector{T}, s :: Selection) where T =
  viewmol(i, x, s.natomspermol)


text:
@testset "Accelerated Gradient Descent" begin
    f(x) = x[1]^4
    function g!(storage, x)
        storage[1] = 4 * x[1]^3
        return
    end

    initial_x = [1.0]
    options = Optim.Options(show_trace = debug_printing, allow_f_increases=true)
    results = Optim.optimize(f, g!, initial_x, AcceleratedGradientDescent(), options)
    @test norm(Optim.minimum(results)) < 1e-6
    @test summary(results) == "Accelerated Gradient Descent"

    # TODO: Check why skip problems fail
    skip = ("Trigonometric", "Large Polynomial", "Parabola", "Paraboloid Random Matrix",
            "Paraboloid Diagonal", "Extended Rosenbrock", "Penalty Function I", "Beale",
            "Extended Powell",
             )
    run_optim_tests(AcceleratedGradientDescent();
                    skip = skip,
                    convergence_exceptions = (("Rosenbrock", 1),("Rosenbrock", 2)),
                    minimum_exceptions = (("Rosenbrock", 2),),
                    minimizer_exceptions = (("Rosenbrock", 2),),
                    iteration_exceptions = (("Powell", 1100),
                                            ("Rosenbrock", 10000),
                                            ("Polynomial", 1500),
                                            ("Fletcher-Powell", 10000),
                                            ("Extended Powell", 8000)),
                    f_increase_exceptions = ("Hosaki", "Polynomial", "Powell", "Himmelblau",
                                             "Extended Powell", "Fletcher-Powell",
                                             "Quadratic Diagonal", "Rosenbrock"),
                    show_name=debug_printing)#,
                    #show_res = debug_printing)
end


text:
using Base.Iterators: repeated, partition

for i in partition(1:60000,32)
    @show typeof(i)
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
struct BVHAccelerator{T,O,B<:BVH{T,O},strategy} <: Accelerator{T,O}
    bvh::B
end
BVHAccelerator(bvh::B,strategy::Symbol) where {T,O,B<:BVH{T,O}} = BVHAccelerator{T,O,B,strategy}(bvh)

function Base.show(io::IO, ::MIME"text/plain", acc::Acc) where {Acc<:BVHAccelerator}
    write(io, "$Acc\n- ")
    show(io, acc.bvh)
end

# * Recursive traversal

const RecursiveBVHAccelerator{T,O,B} = BVHAccelerator{T,O,B,:recursive}

function Base.intersect(ray::Ray{T}, acc::RecursiveBVHAccelerator{T}) where T
    tree = acc.bvh.tree
    node = acc.bvh.root
    intersect(ray, node, tree)
end

function Base.intersect!(intersection::Intersection{T}, acc::RecursiveBVHAccelerator{T}) where T
    tree = acc.bvh.tree
    node = acc.bvh.root
    intersect!(intersection, node, tree)
end

# * Stack traversal

const StackBVHAccelerator{T,O,B} = BVHAccelerator{T,O,B,:stack}

function Base.intersect(ray::Ray{T}, acc::StackBVHAccelerator{T}) where T
    tree = acc.bvh.tree
    node = acc.bvh.root
    intersect(ray, node, tree)
end

function Base.intersect!(intersection::Intersection{T}, acc::StackBVHAccelerator{T}) where T
    tree = acc.bvh.tree
    stack = nodetype(tree)[acc.bvh.root]
    while !isempty(stack)
        node = pop!(stack)
        if is_inner(node)
            if intersect(intersection.ray, node.bbox)
                push!(stack, tree[node.left])
                push!(stack, tree[node.right])
            end
        else
            intersect!(intersection, node)
        end
    end
end

export BVHAccelerator, RecursiveBVHAccelerator, StackBVHAccelerator


text:
@testset "relation.jl" begin
    @testset "sgn" begin
        for i in 0:10
            @test sgn(i) == 1 - nsgn(i) == Int(i != 0)
        end
    end
    @testset "ge" begin
        for i in 0:10, j in rand(0:10, 10)
            @test ge(i, j) == Int(i >= j)
        end
    end
    @testset "ifelse" begin
        for _ in 1:10, b in 0:5
            x, y = rand(0:10, 2)
            @test ifel(b, x, y) == ifelse(b != 0, x, y)
        end
    end
end


text:
using Gadfly

set_default_plot_size(8inch, 3inch)

hstack(
    plot(x=[1,2,3,1],y=[-1,2,3,-4],color=["a","a","a","b"],Geom.bar),
    plot(y=[1,2,3,1],x=[-1,2,3,-4],color=["a","a","a","b"],Geom.bar(orientation=:horizontal)) )


text:
using GeoStatsBase
using DelimitedFiles
using Distances
using Distributions
using Tables, DataFrames
using DensityRatioEstimation
using CategoricalArrays
using LinearAlgebra
using Plots, VisualRegressionTests
using Test, Pkg, Random

using MLJ: @load

# workaround GR warnings
ENV["GKSwstype"] = "100"

# environment settings
islinux = Sys.islinux()
istravis = "TRAVIS" ∈ keys(ENV)
isappveyor = "APPVEYOR" ∈ keys(ENV)
isCI = istravis || isappveyor
visualtests = !isCI || (istravis && islinux)
if !isCI
  Pkg.add("Gtk")
  using Gtk
end
datadir = joinpath(@__DIR__,"data")

# dummy variables for testing
include("dummy.jl")

# list of tests
testfiles = [
  "distances.jl",
  "distributions.jl",
  "data.jl",
  "domains.jl",
  "views.jl",
  "georef.jl",
  "partitioning.jl",
  "weighting.jl",
  "sampling.jl",
  "geoops.jl",
  "paths.jl",
  "trends.jl",
  "geometries.jl",
  "neighborhoods.jl",
  "neighborsearch.jl",
  "learning.jl",
  "mappers.jl",
  "problems.jl",
  "solvers.jl",
  "solutions.jl",
  "errors.jl",
  "statistics.jl",
  "plotrecipes.jl",
  "utils.jl",
  "macros.jl"
]

@testset "GeoStatsBase.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ fabfbc71-7626-4a85-808c-8693d18ddef3
include("..\\..\\jl\\commasplit.jl")

# ╔═╡ 80d3dd1a-1b92-4042-849b-adabbea5a435
include("..\\..\\jl\\arrowsplit.jl")

# ╔═╡ 523ab597-6ce0-40fa-a198-5ff73cabc814
include("..\\..\\jl\\str2int.jl")

# ╔═╡ 610784be-9d2b-47c2-b8b3-623cb0aa9f42
include("..\\..\\jl\\flatten.jl")

# ╔═╡ 87c8c810-657e-11ec-2114-4b4ad4af6ec8
md"""
# Day 5: Hydrothermal Venture

**Reference:** [Day 5 - Advent of Code 2021](https://adventofcode.com/2021/day/5)

**Source:** [ljk233/AdventOfCode2021/days/05/d5.jl](https://github.com/ljk233/AdventOfCode2021/blob/main/days/05/d5.jl)
"""

# ╔═╡ b86823bf-1386-449f-8940-17aa45a0d242
md"""
## Imports
"""

# ╔═╡ e05f10b5-f57e-42ca-8ee6-5df06c0baa59
md"""
## Functions
"""

# ╔═╡ 46dac4eb-b5d9-4fa5-badb-9d69f223c883
"""
	V(line: String) -> NamedTuple(NamedTuble, NamedTuple)
"""
V(lines) = map(lines) do line
	r = line |> arrowsplit .|> commasplit |> flatten .|> str2int
	return (a=(x=r[1]+1, y=r[2]+1), b=(x=r[3]+1, y=r[4]+1))
end

# ╔═╡ a5a42a42-4601-408e-920a-72cf7c818ac2
"""
	mark!(grid, a, b, [inc45=false])
"""
function mark!(grid, a, b, inc45=false)
	if inc45 && abs((a.y-b.y)/(a.x-b.x)) == 1
		dx = 0
		if a.x > b.x
			a, b = b, a
		end
		for y ∈ a.y:cmp(b.y, a.y):b.y
			grid[a.x+dx, y] += 1
			dx+=1
		end
	elseif a.x == b.x
		for y ∈ a.y:cmp(b.y, a.y):b.y
			grid[a.x, y] += 1
		end
	elseif a.y == b.y
		for x ∈ a.x:cmp(b.x, a.x):b.x
			grid[x, a.y] += 1
		end
	end
end

# ╔═╡ 4c202df3-62b5-4446-b3e1-0afae48210b8
md"""
## Part One

To avoid the most dangerous areas, you need to determine the number of points where at least two lines overlap.
Consider only horizontal and vertical lines.

**At how many points do at least two lines overlap?**
"""

# ╔═╡ 1a9694c2-923b-4f98-83f6-2b2bf3ff55ab
"""
	main1(f: String) -> Integer
"""
function main1(f)
	lines = readlines(f)
	xmax = maximum(v -> max(v.a.x, v.b.x), V(lines))
	ymax = maximum(v -> max(v.a.x, v.b.x), V(lines))
	grid = zeros(Integer, (xmax+1, ymax+1))
	for (a, b) ∈ V(lines)
		mark!(grid, a, b)
	end
	return count(>(1), grid)
end

# ╔═╡ c2a8fd18-c49e-48e0-8780-0591a3c59fa6
md"""
**Test #1.**
Expected result = 5.
"""

# ╔═╡ a5ac2615-cae0-4a16-901d-687fd46f3145
main1("test.in")

# ╔═╡ a957d84e-b1a4-4a95-8539-fcaf07f1f1e5
md"""
**Solution.**
Confirmed result = 8350.
"""

# ╔═╡ d6b93da4-1ef1-48d6-83fe-8424acbcfaf2
main1("data.in")

# ╔═╡ cede829f-55f0-4f4d-96b5-a380acd16394
md"""
## Part Two

You still need to determine the number of points where at least two lines overlap.

**Consider all of the lines.**
**At how many points do at least two lines overlap?**
"""

# ╔═╡ 3bc8ca5f-ba0a-4943-a266-59e2079365cc
"""
	main2(f: String) -> Integer
"""
function main2(f)
	lines = readlines(f)
	xmax = maximum(v -> max(v.a.x, v.b.x), V(lines))
	ymax = maximum(v -> max(v.a.x, v.b.x), V(lines))
	grid = zeros(Integer, (xmax+1, ymax+1))
	for (a, b) ∈ V(lines)
		mark!(grid, a, b, true)
	end
	return count(>(1), grid)
end

# ╔═╡ 2d234229-dbea-40cc-be6e-0267a04bd1b0
md"""
**Test #2.**
Expected result = 12.
"""

# ╔═╡ e099ec40-a9fe-4446-9db7-9619cb3adc5d
main2("test.in")

# ╔═╡ a9688396-74a2-4989-95bf-381b7904f521
md"""
**Solution.**
Confirmed result = 19374.
"""

# ╔═╡ 488943e4-d1ba-4818-b8d9-1f886dc7529f
main2("data.in")

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╟─87c8c810-657e-11ec-2114-4b4ad4af6ec8
# ╟─b86823bf-1386-449f-8940-17aa45a0d242
# ╠═fabfbc71-7626-4a85-808c-8693d18ddef3
# ╠═80d3dd1a-1b92-4042-849b-adabbea5a435
# ╠═523ab597-6ce0-40fa-a198-5ff73cabc814
# ╠═610784be-9d2b-47c2-b8b3-623cb0aa9f42
# ╟─e05f10b5-f57e-42ca-8ee6-5df06c0baa59
# ╠═46dac4eb-b5d9-4fa5-badb-9d69f223c883
# ╠═a5a42a42-4601-408e-920a-72cf7c818ac2
# ╟─4c202df3-62b5-4446-b3e1-0afae48210b8
# ╠═1a9694c2-923b-4f98-83f6-2b2bf3ff55ab
# ╟─c2a8fd18-c49e-48e0-8780-0591a3c59fa6
# ╠═a5ac2615-cae0-4a16-901d-687fd46f3145
# ╟─a957d84e-b1a4-4a95-8539-fcaf07f1f1e5
# ╠═d6b93da4-1ef1-48d6-83fe-8424acbcfaf2
# ╟─cede829f-55f0-4f4d-96b5-a380acd16394
# ╠═3bc8ca5f-ba0a-4943-a266-59e2079365cc
# ╟─2d234229-dbea-40cc-be6e-0267a04bd1b0
# ╠═e099ec40-a9fe-4446-9db7-9619cb3adc5d
# ╟─a9688396-74a2-4989-95bf-381b7904f521
# ╠═488943e4-d1ba-4818-b8d9-1f886dc7529f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module PT1

ODE-model of PT1 block (T*dx/dt + x = 0).
"""
module PT1

import ModiaMath

mutable struct Model <: ModiaMath.AbstractSimulationModel
    simulationState::ModiaMath.SimulationState

    # Parameters
    T::Float64

    # Initial conditions
    x0::Float64

    function Model(;T=1.0, x0=1.0)
        @assert(T > 0.0)
        simulationState = ModiaMath.SimulationState("PT1", getModelResidues!, [x0];
                                                    structureOfDAE = ModiaMath.DAE_ExplicitDerivatives)
        new(simulationState, T, x0)
    end
end

function getModelResidues!(m::Model, t::Float64, _x::Vector{Float64}, _derx::Vector{Float64}, _r::Vector{Float64}, _w::Vector{Float64})
    x     = _x[1]
    derx  = -x / m.T
    _r[1] = derx - _derx[1]
    return nothing
end

end


text:
module IncrementalInference

import Base: convert
import HDF5: root
# import KernelDensityEstimate: root
import Distributions: sample
# import KernelDensityEstimate: sample
import KernelDensityEstimate: kde!
# import Graphs: plot

using
  Graphs,
  # GraphViz,
  NLsolve,
  Optim,
  Distributions,
  KernelDensityEstimate,
  HDF5,
  JLD,
  ProgressMeter,
  DocStringExtensions,
  Compat

export
  # pass through from Graphs.jl
  # plot,

  # added methods to functions from KernelDensityEstimate
  kde!,

  # pass through functions commonly used lower down
  Npoints,
  Ndim,
  getBW,

  # data layer variables
  dlapi,
  localapi,
  showcurrentdlapi,
  # and callback setting function
  setdatalayerAPI!,
  DataLayerAPI,

  # general types for softtyping of variable nodes
  InferenceVariable,
  ContinuousScalar,
  ContinuousMultivariate,

  # using either dictionary or cloudgraphs
  VariableNodeData,
  PackedVariableNodeData,
  VNDencoder,
  VNDdecoder,
  FNDencode,
  FNDdecode,
  FunctionNodeData,
  PackedFunctionNodeData,
  FactorGraph,
  addNode!,
  addFactor!,
  resetData!,
  getVert,
  getData,
  setData!,
  getVarNode,
  getVal,
  setVal!,
  getBWVal,
  setBW!,
  setValKDE!,
  isInitialized,
  ensureAllInitialized!,
  updateFullVert!,
  getOutNeighbors,
  BayesTree,
  EasyMessage,
  NBPMessage,
  ExploreTreeType,
  emptyFactorGraph,
  subgraphShortestPath,
  subgraphFromVerts,
  getEliminationOrder,
  buildBayesNet!,
  emptyBayesTree,
  buildTree!,
  prepBatchTree!,
  wipeBuildNewTree!,
  whichCliq,
  getKDE,
  getVertKDE,
  initializeNode!,

  #functors need
  getSample,

  #Visualization
  # writeGraphPdf,
  ls,
  lsf,
  ls2,
  getfnctype,
  drawCopyFG,
  # investigateMultidimKDE,
  # drawHorDens,
  # drawHorBeliefsList,

  # Tree stuff
  # spyCliqMat,
  evalPotential,
  evalFactor2,

  # dev
  GenericWrapParam,

  # solve inference
  inferOverTree!,
  inferOverTreeR!,
    #development interface
    upMsgPassingRecursive,

  # Inference types
  InferenceType,
  PackedInferenceType,
  Singleton,
  Pairwise,

  # functor abstracts
  FunctorInferenceType,
  FunctorPairwise,
  FunctorSingleton,
  # FunctorPartialSingleton,
  FunctorPairwiseNH,
  FunctorSingletonNH,

  # Solving utils
  findRelatedFromPotential,
  shuffleXAltD!,
  numericRoot,
  numericRootGenericRandomized,
  numericRootGenericRandomizedFnc,
  numericRootGenericRandomizedFnc!,
  FastGenericRoot,
  FastRootGenericWrapParam,

  # user functions
  proposalbeliefs,
  predictbelief,

  # generic marginal used during elimitation game
  GenericMarginal,
  PackedGenericMarginal,

  uppA,
  convert, # for protobuf stuff
  compare,

  # factor graph operating system utils (fgos)
  convert2packedfunctionnode,
  encodefg,
  decodefg,
  savejld,
  loadjld,
  landmarks,

  # For 1D example, should be refactored and renamed
  Odo,
  odoAdd,
  PackedOdo,
  Obsv2,
  PackedObsv2,
  Ranged,
  PackedRanged,

  # dev exports
  addGraphsVert!,
  makeAddEdge!,

  # define evalPotential functions outside IIF
  registerCallback!,

  # development
  shuffleXAltD,
  reshapeVec2Mat

  # analysis and some plotting
  # plotKDEofnc,
  # plotKDEresiduals,
  # plotMCMC,
  # plotUpMsgsAtCliq,
  # plotPriorsAtCliq,
  # investigateMultidimKDE,
  # draw,
  # whosWith,
  # drawUpMsgAtCliq,
  # dwnMsgsAtCliq,
  # drawPose2DMC!,
  # mcmcPose2D!,
  # # drawUpMCMCPose2D!,
  # # drawDwnMCMCPose2D!,
  # drawLbl,
  # predCurrFactorBeliefs,
  # drawHorDens,
  # drawHorBeliefsList,
  # drawFactorBeliefs,
  # localProduct,
  # plotLocalProduct,
  # saveplot,
  # animateVertexBelief,
  # getColorsByLength

const VoidUnion{T} = Union{Void, T}

include("FactorGraphTypes.jl")
include("DataLayerAPI.jl")
include("FactorGraph01.jl")
include("JunctionTree.jl")
include("GraphConstraintTypes.jl")
include("SolverUtilities.jl")
include("TreePotentials01.jl")
include("ApproxConv.jl")
include("SolveTree01.jl")
#include("SolverVisualization.jl")
include("FGOSUtils.jl")

include("deprecated.jl")

# function plot(fg::FactorGraph)
#   Graphs.plot(fg.g)
# end

end


text:
using Distributions
export levitt_adapter, bayesian_adapter, constant_adapter, delta, estimate,
  update!

abstract type Adapter end

"""
    update!(adapter,response,correct)

Updates any internal state for the adapter when the listener responds with
`response` and the correct response is `correct`. Usually not called directly,
but instead called within `response`, when the adapter is passed as the first
argument. May take a while to run.
"""
function update!(adapter,response,correct)
end

"""
    estimate(adapter)

Returns the mean and error of the adapters threshold estimate. May take some
time to run.
"""
function estimate(adapter)
end

"""
    delta(adapter)

Returns the next delta that should be tested to help estimate the threshold.
"""
function delta(adapter)
end

"""
    response([fn],adapter,[key1] => ["resp1"],...;correct=[resp],
             [show_feedback=true],
             [feedback=Dict(true=>"Correct",false=>"Wrong!")]
             keys...)

Record a response in a n-alternative forced choice task and update
an adapter.

The first response recieved is interpreted as the actual response. Subsequent
responses will be recorded, without a delta or correct value set, and appending
"late_" to the specified response string.

# Function Callback

Optionally, upon participant response, `fn` receives two arguments: the
provided response, and the correct response.

# Keyword Arguments

- `correct` the response string corresponding to the correct response
- `show_feedback` (default = true): whether to show feedback to the
  participant after they respond.
- `feedback` the text to display to a participant when they are correct
  (for the true key, defaults to "Correct!") or incorrect (for the false key,
  defaults to "Wrong!").

Any additional keyword arguments are added as column values when the
response is recorded.

"""
function response(track::Adapter,resp::Pair...;keys...)
  response((x,y) -> nothing,track,resp...;keys...)
end

function response(callback::Function,adapter::Adapter,responses::Pair...;
                  correct=nothing,show_feedback=true,
                  feedback=Dict(true => "Correct!",false => "Wrong!"),
                  keys...)
  addcolumn(:correct)
  addcolumn(:delta)

  if correct ∉ map(x -> x[2],responses)
    error("The value of `correct` must be "*
          join(map(x -> x[2],responses),", "," or "))
  end
  feedback_visuals = map((pair) -> pair[1] => visual(pair[2]),feedback)

  let responded = false
    begin (event) ->
      for (key,resp) in responses
        if iskeydown(event,key)
          if !responded
            responded = true
            record(resp;correct=correct,delta=delta(adapter),keys...)

            update!(adapter,resp,correct)
            callback(resp,correct)
            if show_feedback
              display(feedback_visuals[resp==correct])
            end

          else
            record("late_"*resp;keys...)
          end
        end
      end
    end
  end
end

mutable struct ConstantStimulus
  stimuli::Vector{Float64}
  index::Int
  correct::Vector{Int}
  response::Vector{Int}
end

"""
    constant_adapter(stimuli)

An adapter that can be used to implement the method of constant stimuli: the
specified sequence of stimulus deltas is presented in order to participants.

Strictly speaking, this is not an adaptive tracking procedure. However,
it can be convienient to have the same programming interface for this method as
for adaptive methods. In this way you can easily select between the method of
constant stimuli or some kind of adaptive procedure.
"""
constant_adapter(stimuli) = ConstantStimulus(stimuli,1,Int[],Int[])
delta(adapter::ConstantStimulus) = adapter.stimuli[adapter.index]
function update!(adapter::ConstantStimulus,response,correct)
  adapte.index += 1
  nothing
  # push!(adapter.correct,correct)
  # push!(adapter.response,response)
end
function estimate(adapter::ConstantStimulus)
  error("The `estimate` method for `ConstantStimulus` is not implemented!")
  # TODO: fit a curve to the responses to estimate the threshold???
end

"""
    levitt_adapter([first_delta=0.1],[down=3],[up=1],
                   [big_reverse=3],[big=0.01],[little=0.005],
                   [min_reversals=7],[min_delta=-Inf],[max_delta=Inf],
                   [mult=false])

An adapter that finds a threshold according to a non-parametric statistical
procedure. This approach makes fewer explicit assumptions than
[`bayesian_adapter`](@ref) but may be slower to converge to a threshold.

This finds a threshold by moving the delta down after three correct responses
and up after one incorrect response (these default up and down counts can be
changed). This is the same approach described in Levitt (1971).

# Keyword Arguments

- `first_delta`: the delta that the first trial should present.
- `up`: how many incorrect responses in a row must occur for the delta
  to move up
- `down`: how many correct responses in a row must occur for the delta
  to move down.
- `big`: the amount delta changes by (up or down) at first
- `big_reverse`: how many reveresals (up to down or down to up) must
  occur before `little` is used instead of `big`
- `little`: the amount delta changes by (up or down) after
  `big_reverse` reversals.
- `min_reversals`: the smallest number of reversals that can
  be used to estimate a threshold.
- `min_delta`: the smallest delta allowed
- `max_delta`: the largest delta allowed
- `mult`: whether the delta change should be additive (false) or
   multiplicative (true).
"""
function levitt_adapter(;first_delta=0.1,down=3,up=1,big_reverse=3,
                        big=0.01,little=0.005,min_reversals=7,
                        drop_reversals=3,min_delta=0,max_delta=Inf,
                        mult=false)
  delta = first_delta
  reversals = Float64[]
  last_direction = 0
  num_correct = 0
  num_incorrect = 0
  operator = (mult ? :mult : :add)
  adapter = Levitt{operator}(delta,up,down,big_reverse,big,little,
                             min_reversals,drop_reversals,min_delta,max_delta,
                             reversals,last_direction,num_correct,num_incorrect)
end

type Levitt{OP} <: Adapter
  delta::Float64

  up::Int
  down::Int
  big_reverse::Int
  big::Float64
  little::Float64
  min_reversals::Int
  drop_reversals::Int

  min_delta::Float64
  max_delta::Float64

  reversals::Array{Float64}
  last_direction::Int
  num_correct::Int
  num_incorrect::Int
end
delta(a::Levitt) = a.delta

function update_reversals(adapter::Levitt,direction::Int)
  if (adapter.last_direction != 0 && direction != adapter.last_direction)
    push!(adapter.reversals,adapter.delta)
  end
  adapter.last_direction = direction
end

bound(x,min_,max_) = min(max_,max(min_,x))

up(adapter::Levitt{:add}) = +
down(adapter::Levitt{:add}) = -
up(adapter::Levitt{:mult}) = *
down(adapter::Levitt{:mult}) = /

function update_delta(adapter,op)
  if length(adapter.reversals) < adapter.big_reverse
    adapter.delta = bound(op(adapter.delta,adapter.big),
                          adapter.min_delta,adapter.max_delta)
  else
    adapter.delta = bound(op(adapter.delta,adapter.little),
                          adapter.min_delta,adapter.max_delta)
  end
end

function update!(adapter::Levitt,response,correct)
  if correct == response
    adapter.num_correct += 1
    adapter.num_incorrect = 0
    if adapter.num_correct >= adapter.down
      adapter.num_correct = 0
      update_reversals(adapter,-1)
      update_delta(adapter,down(adapter))
    end
  else
    adapter.num_incorrect += 1
    adapter.num_correct = 0
    if adapter.num_incorrect >= adapter.up
      adapter.num_incorrect = 0
      update_reversals(adapter,1)
      update_delta(adapter,up(adapter))
    end
  end
  adapter
end

estimate_helper(a::Levitt{:add},rev) = mean(rev),std(rev)
estimate_helper(a::Levitt{:mult},rev) =
  exp(mean(log(rev))),exp(mean(log(rev)))*(1-exp(-std(log(rev))))

function estimate(adapter::Levitt)
  if length(adapter.reversals) < adapter.min_reversals
    NaN,NaN
  else
    if isodd(length(adapter.reversals))
      estimate_helper(adapter,adapter.reversals[adapter.drop_reversals:end])
    else
      estimate_helper(adapter,adapter.reversals[adapter.drop_reversals+1:end])
    end
  end
end

type ImportanceSampler <: Adapter
  miss::Float64
  threshold::Float64
  thresh_d::UnivariateDistribution
  inv_slope_d::UnivariateDistribution
  n_samples::Int

  theta_prior::UnivariateDistribution
  inv_slope_prior::UnivariateDistribution

  delta::Float64
  min_delta::Float64
  max_delta::Float64
  delta_repeat::Int

  resp::Dict{Float64,Int}
  N::Dict{Float64,Int}
  thresh_samples::Vector{Float64}
  theta_samples::Vector{Float64}
  inv_slope_samples::Vector{Float64}
  thresh_weights::Vector{Float64}

  repeat2_thresh::Float64
  repeat3_thresh::Float64
end

# TODO: prior is working, but improve it to make the
# parameters easy to specify???

"""
    bayesian_adapter(;first_delta=0.1,
                     n_samples=1000,miss=0.01,threshold=0.79,
                     min_delta=0,max_delta=1,
                     min_plausible_delta = 0.0001,
                     max_plausible_delta = 0.2,
                     repeat3_thresh=1.0,repeat2_thresh=0.1,
                     thresh_prior=
                     Truncated(LogNormal(log(min_plausible_delta),
                                         log(max_plausible_delta/
                                             min_plausible_delta/2)),
                               min_delta,max_delta),
                          inv_slope_prior=TruncatedNormal(0,0.25,0,Inf),
                          thresh_d=thresh_prior,
                          inv_slope_d=inv_slope_prior)

An adapter that finds a threshold according to a parametric statistical
model. This makes more explicit assumptions than the [`levitt_adapter`](@ref)
but will normally find a threshold faster.

The psychometric curve is estimated from user responses using a bayesian
approach. After estimation, each new delta is selected in a greedy fashion,
picking the response that best minimizes entropy according to this psychometric
function. This is a modified version of the approach described in Kontsevich &
Tyler 1999. Specifically, the approach here uses importance sampling instead of
sampling parmeters over a deterministic, uniform grid. This should increase
measurement efficiency if the priors are chosen well.

This algorithm assumes the following functional form for the psychometric
response as a function of the stimulus difference ``Δ``.

``
f(Δ) = λ/2 + (1-λ) Φ((Δ - θ)⋅σ/√2)
``

In the above ``Φ`` is the cumulative distribution function of a normal
distribution, ``λ`` is the miss-rate parameter, indicating the rate at which
listeners make a mistake, even when the delta is large and easily heard, ``θ``
is the 50%-correct threshold, and ``σ`` is the psychometric slope.

For stability and robustness, this adapter begins by repeating the same delta
multiple times and only begins quickly changing deltas trial-by-trial when the
ratio of estiamted standard deviation to mean is small. This functionality can
be adjusted using `repeat3_thresh` and `repeat2_thresh`, or, if you do not wish
to have any repeats, both values can be set to Inf.

# Keyword Arugments

- **first_delta**: the delta to start measuring with
- **n_samples** the number of samples to use during importance sampling.
  The algorithm for selecting new deltas is O(n²).
- **miss** the expected rate at which listeners will make mistakes
  even for easy to percieve differences.
- **threshold** the %-response threshold to be estimated
- **min_delta** the smallest possible delta
- **max_delta** the largest possible delta
- **min_plausible_delta** the smallest plausible delta, should be > 0.
  Used to define a reasonable value for thresh_prior and inv_slope_prior.
- **max_plausible_delta** the largest plausible delta, should be < max_delta.
  Used to define a reasonable value for thresh_prior and inv_slope_prior.
- **thresh_prior** the prior probability distribution across thresholds.
  This influence the way the delta is adapted. By default this is defined in
  terms of min_plausible_delta and max_plausible_delta.
- **inv_slope_prior** the prior probability distribution across inverse slopes.
  By default this is defined in terms of min_plausible_delta and
  max_plausible_delta.
- **thresh_d** the distribution over-which to draw samples for the threshold
  during importance sampling. This defaults to thresh_prior
- **inv_slope_d** the distribution over-which to draw samples for the inverse slope
  during importance sampling. This defaults to inv_slope_prior.
- **repeat2_thresh** the ratio of sd / mean for theta must suprass
  to repeat each delta twice.
- **repeat3_thresh** the ratio of sd / mean for theta must surpass
  to repeat each delta thrice.
"""
function bayesian_adapter(;first_delta=0.1,
                          n_samples=1000,miss=0.01,threshold=0.79,
                          min_delta=0,max_delta=1,
                          min_plausible_delta = 0.0001,
                          max_plausible_delta = 0.2,
                          repeat3_thresh=1.0,repeat2_thresh=0.1,
                          thresh_prior=
                          Truncated(LogNormal(log(min_plausible_delta),
                                              log(max_plausible_delta/
                                                  min_plausible_delta/2)),
                                    min_delta,max_delta),
                          inv_slope_prior=
                          TruncatedNormal(0,2max_plausible_delta,0,Inf),
                          thresh_d=thresh_prior,
                          inv_slope_d=inv_slope_prior)
  delta = first_delta
  delta_repeat = 0
  resp = Dict{Float64,Int}()
  N = Dict{Float64,Int}()
  samples = ones(0)
  weights = ones(0)
  ImportanceSampler(miss,threshold,thresh_d,inv_slope_d,n_samples,
                    thresh_prior,inv_slope_prior,
                    delta,min_delta,max_delta,delta_repeat,
                    resp,N,samples,samples,samples,weights,
                    repeat3_thresh,repeat2_thresh)
end

const sqrt_2 = sqrt(2)
function sample(a::ImportanceSampler)
  if length(a.thresh_samples) == 0
    theta = rand(a.thresh_d,a.n_samples)
    inv_slope = rand(a.inv_slope_d,a.n_samples)
    q_log_prob = logpdf(a.thresh_d,theta) .+ logpdf(a.inv_slope_d,inv_slope)

    udeltas = collect(keys(a.resp))
    x = exp((udeltas' .- theta)./inv_slope)/sqrt_2
    p = (a.miss/2) + (1-a.miss)*cdf(Normal(),x)
    log_prob = sum(enumerate(udeltas)) do x
      i,delta = x
      logpdf.(Binomial.(a.N[delta],p[:,i]),a.resp[delta])
    end + logpdf(a.theta_prior,theta) + logpdf(a.inv_slope_prior,inv_slope)

    lweights = log_prob - q_log_prob
    weights = exp.(lweights - maximum(lweights))
    weights = weights ./ sum(weights)
    ess = 1 / sum(weights.^2)

    # if ess / length(weights) < 0.1
    #   warn("Effective sampling size of importance samples is low ",
    #        "($(round(ess,1))). Consider adjusting the `thresh_d` and `slope_d`.")
    #   record("poor_adaptor_ess",value=ess)
    # end

    dprime = invlogcdf(Normal(),log((a.threshold-a.miss/2)/(1-a.miss)))
    thresh_off = -log(dprime).*inv_slope
    a.theta_samples = theta
    a.inv_slope_samples = inv_slope
    a.thresh_samples = theta + thresh_off
    a.thresh_weights = weights
  end
  a.thresh_samples, a.thresh_weights
end

function update!(adapter::ImportanceSampler,response,correct)
  adapter.thresh_samples = ones(0)
  adapter.resp[adapter.delta] = get(adapter.resp,adapter.delta,0) +
    (response==correct)
  adapter.N[adapter.delta] = get(adapter.N,adapter.delta,0) + 1

  thresh,weights = sample(adapter)
  μ,σ = estimate(adapter)

  sd_ratio = σ / μ
  if sd_ratio > adapter.repeat3_thresh
    repeats = 2
  elseif sd_ratio > adapter.repeat2_thresh
    repeats = 1
  else
    repeats = 0
  end

  if repeats > adapter.delta_repeat
    adapter.delta_repeat += 1
  else
    adapter.delta_repeat = 0

    # TODO: sample these differences to make sure this is a reasonably
    # sized array
    # TODO: debug selection

    entropies = -sum(1:adapter.n_samples) do i
      diff = adapter.thresh_samples[i] .- adapter.theta_samples
      dprime = exp(diff./adapter.inv_slope_samples)/sqrt_2
      p = ((adapter.miss/2) + (1-adapter.miss)*cdf(Normal(),dprime))
      (p .* log(p) + (1-p) .* log(1-p)) .* weights
    end

    _,i = findmax(entropies)

    adapter.delta = adapter.thresh_samples[i]
  end

  adapter
end

function estimate(adapter::ImportanceSampler)
  thresh,weights = sample(adapter)
  μ = sum(thresh .* weights)
  M = sum(weights .> 0)
  σ = sqrt.(sum(((μ .- thresh) .* weights).^2) ./
           ((M-1)/M)*sum(weights))

  μ,σ
end

delta(adapter::ImportanceSampler) = bound(adapter.delta,adapter.min_delta,
adapter.max_delta)


text:
# Code for toy example involving Bernoulli trials
module bernoulli_example
using PyPlot, Distributions
draw_now() = (pause(0.001); get_current_fig_manager()[:window][:raise_]())
logsumexp(x) = (m = maximum(x); m == -Inf ? -Inf : log(sum(exp(x-m))) + m)

# H0: theta = 1/2
# H1: theta \neq 1/2, theta|H1 ~ Beta(1,1)

# settings
epsilon = 0.02 # "precision" of theta
alpha = 1/(2*epsilon^2) # choose alpha using relative entropy approach, for theta near 1/2
theta0 = 0.56 # true value
nns = 25 # 1+6*4 = number of n's to use
nreps = 1000 # number of times to run the simulation
ns = int(logspace(0,6,nns)) # n's to use

# initialize
p0xs = zeros(nreps,nns)
r0xs = zeros(nreps,nns)
e0xs = zeros(nreps,nns)
maxn = maximum(ns)

for rep = 1:nreps
    y = (rand(maxn).<theta0) # data
    for (i,n) in enumerate(ns)
        s = sum(y[1:n])

        # standard posterior
        B10 = exp(n*log(2) + lbeta(1+s, 1+n-s)) # p(x|H1)/p(x|H0)
        p0x = 1/(1+B10)  # p(H0|x)

        # robust posterior - approximate
        alpha_n = 1/(1/alpha + 1/n)
        R10 = exp(alpha_n*log(2) + lbeta(1+alpha_n*s/n, 1+alpha_n*(1-s/n)))
        r0x = 1/(1+R10)  # p(H0|E)

        # robust posterior - exact
        f(p,q) = (r=p.*log(p./q); if length(p)>1; r[p.==0]=0; elseif p==0; r[:]=0; end; r)
        D(p,q) = f(p,q) + f(1-p,1-q) # relative entropy for Bernoulli
        t = [0:n]
        lpt0 = logpdf(Binomial(n,1/2),t) # Binomial(n,1/2)
        lpt1 = -log(n+1)*ones(n+1) # BetaBinomial(n,1,1) = Uniform{0,1,...,n}
        lpE0 = logsumexp(-alpha*D(s/n,t/n) + lpt0)  # log(p(E|H0))
        lpE1 = logsumexp(-alpha*D(s/n,t/n) + lpt1)  # log(p(E|H1))
        E10 = exp(lpE1 - lpE0)
        e0x = 1/(1+E10)  # p(H0|E)

        # record values
        p0xs[rep,i] = p0x
        r0xs[rep,i] = r0x
        e0xs[rep,i] = e0x

        #println("n = $n")
        #@printf("p(H0|x)  std=%.4f  approx=%.5f  exact=%.5f\n",p0x,r0x,e0x)
    end
end

# plot
figure(1,figsize=(10,3.2)); clf(); hold(true); grid(true)
subplots_adjust(bottom = 0.2)
semilogx(ns,mean(p0xs,1)[:],label="standard posterior",linewidth=2,"gs-",markersize=5)
semilogx(ns,mean(e0xs,1)[:],label="exact c-posterior",linewidth=2,"r*-",markersize=10)
semilogx(ns,mean(r0xs,1)[:],label="approx c-posterior",linewidth=1,"bo-",markersize=5)
xlabel("\$n \\,\\,\\mathrm{(sample\\,size)}\$",fontsize=16)
ylabel("\$\\Pi(\$H\$_0|\\mathrm{data})\$",fontsize=15)
ylim(0,1)
#legend(numpoints=1,loc="lower left",fontsize=14)
draw_now()
savefig("pH0-$theta0.png",dpi=150)

end

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
using GLPK, Base.Test, MathOptInterface, MathOptInterface.Test, MathOptInterfaceGLPK

const MOIT = MathOptInterface.Test
const MOIGLP = MathOptInterfaceGLPK

@testset "MathOptInterfaceGLPK" begin
    @testset "Linear tests" begin
        linconfig_nocertificate = MOIT.TestConfig(infeas_certificates=false)
        solver = GLPKOptimizerLP()
        MOIT.contlineartest(solver, linconfig_nocertificate, ["linear10"])

        linconfig_nocertificate_noduals = MOIT.TestConfig(duals=false,infeas_certificates=false)
        solver_mip = GLPKOptimizerMIP()
        MOIT.contlineartest(solver_mip, linconfig_nocertificate_noduals, ["linear10","linear8b","linear8c"])

        # Intervals
        linconfig_noquery = MOIT.TestConfig(query=false)
        MOIT.linear10test(solver, linconfig_noquery)
    end

    @testset "Linear Conic tests" begin
        linconfig_nocertificate = MOIT.TestConfig(infeas_certificates=false)
        solver = GLPKOptimizerLP()
        MOIT.lintest(solver, linconfig_nocertificate)

        linconfig_nocertificate_noduals = MOIT.TestConfig(duals=false,infeas_certificates=false)
        solver_mip = GLPKOptimizerMIP()
        MOIT.lintest(solver_mip, linconfig_nocertificate_noduals)
    end

    @testset "Integer Linear tests" begin
        intconfig = MOIT.TestConfig()
        solver_mip = GLPKOptimizerMIP()
        MOIT.intlineartest(solver_mip, intconfig, ["int2","int1"])
    end

    @testset "ModelLike tests - MIP" begin
        intconfig = MOIT.TestConfig()
        solver = GLPKOptimizerMIP()
        MOIT.validtest(solver)
        MOIT.emptytest(solver)
        solver2 = GLPKOptimizerMIP()
        MOIT.copytest(solver,solver2)
    end
    @testset "ModelLike tests - LP" begin
        intconfig = MOIT.TestConfig()
        solver = GLPKOptimizerLP()
        MOIT.validtest(solver)
        MOIT.emptytest(solver)
        solver2 = GLPKOptimizerLP()
        MOIT.copytest(solver,solver2)
    end
end
;


text:
using Test, Distributions, ForwardDiff

isnan_type(::Type{T}, v) where {T} = isnan(v) && v isa T

@testset "Normal" begin
    @test isa(convert(Normal{Float64}, Float16(0), Float16(1)),
              Normal{Float64})
    @test Inf === logpdf(Normal(0, 0), 0)
    @test -Inf === logpdf(Normal(), Inf)
    @test iszero(logcdf(Normal(0, 0), 0))
    @test iszero(logcdf(Normal(), Inf))
    @test logdiffcdf(Normal(), Float32(5), Float32(3)) ≈ -6.607938594596893 rtol=1e-12
    @test logdiffcdf(Normal(), Float32(5), Float64(3)) ≈ -6.607938594596893 rtol=1e-12
    @test logdiffcdf(Normal(), Float64(5), Float64(3)) ≈ -6.607938594596893 rtol=1e-12
    let d = Normal(Float64(0), Float64(1)), x = Float64(-60), y = Float64(-60.001)
        float_res = logdiffcdf(d, x, y)
        big_float_res = log(cdf(d, BigFloat(x, 100)) - cdf(d, BigFloat(y, 100)))
        @test float_res ≈ big_float_res
    end
    @test_throws ArgumentError logdiffcdf(Normal(), 1.0, 2.0)
    @test -Inf === logccdf(Normal(0, 0), 0)
    @test iszero(logccdf(Normal(eps(), 0), 0))
    @test -Inf === quantile(Normal(), 0)
    @test iszero(quantile(Normal(), 0.5))
    @test Inf === quantile(Normal(), 1)
    @test -Inf === quantile(Normal(0, 0), 0)
    @test iszero(quantile(Normal(0, 0), 0.75))
    @test Inf === quantile(Normal(0, 0), 1)
    @test -Inf === quantile(Normal(0.25, 0), 0)
    @test 0.25 == quantile(Normal(0.25, 0), 0.95)
    @test Inf === quantile(Normal(0.25, 0), 1)
    @test Inf === cquantile(Normal(), 0)
    @test iszero(cquantile(Normal(), 0.5))
    @test -Inf === cquantile(Normal(), 1)
    @test Inf === cquantile(Normal(0, 0), 0)
    @test iszero(cquantile(Normal(0, 0), 0.75))
    @test -Inf === cquantile(Normal(0, 0), 1)
    @test Inf === cquantile(Normal(0.25, 0), 0)
    @test 0.25 == cquantile(Normal(0.25, 0), 0.95)
    @test -Inf === cquantile(Normal(0.25, 0), 1)
    @test -Inf === invlogcdf(Normal(), -Inf)
    @test isnan_type(Float64, invlogcdf(Normal(), NaN))
    @test Inf === invlogccdf(Normal(), -Inf)
    @test isnan_type(Float64, invlogccdf(Normal(), NaN))
    # test for #996 being fixed
    let d = Normal(0, 1), x = 1.0, ∂x = 2.0
        @inferred cdf(d, ForwardDiff.Dual(x, ∂x)) ≈ ForwardDiff.Dual(cdf(d, x), ∂x * pdf(d, x))
    end
end

@testset "Normal logpdf & pdf type inference" begin
    @test @inferred(pdf(Normal(0.0, 0.0), 0.0))           === Inf
    @test @inferred(pdf(Normal(0.0, 0.0), -1.0))          === 0.0
    @test @inferred(pdf(Normal(0.0, 0.0), 0.0f0))         === Inf
    @test isnan_type(Float64, @inferred(pdf(Normal(0.0, 0.0), NaN)))
    @test @inferred(pdf(Normal(0.0f0, 0.0f0), 0.0))       === Inf
    @test @inferred(pdf(Normal(0.0f0, 0.0f0), 0.0f0))     === Inf32
    @test isnan_type(Float64, @inferred(pdf(Normal(0.0, 0.0), NaN)))
    @test isnan_type(Float64, @inferred(pdf(Normal(NaN, 0.0), 0.0f0)))
    @test isnan_type(Float32, @inferred(pdf(Normal(NaN32, 0.0f0), 0.0f0)))
    @test @inferred(pdf(Normal(0 // 1, 0 // 1), 0 // 1))  === Inf
    @test isnan_type(Float64, @inferred(pdf(Normal(0 // 1, 0 // 1), NaN)))
    @test @inferred(pdf(Normal(0.0, 0.0), BigInt(1)))     == big(0.0)
    @test @inferred(pdf(Normal(0.0, 0.0), BigFloat(1)))   == big(0.0)
    @test isnan_type(BigFloat, @inferred(pdf(Normal(0.0, 0.0), BigFloat(NaN))))

    @test @inferred(logpdf(Normal(0.0, 0.0), 0.0))           === Inf
    @test @inferred(logpdf(Normal(0.0, 0.0), -1.0))          === -Inf
    @test @inferred(logpdf(Normal(0.0, 0.0), 0.0f0))         === Inf
    @test isnan_type(Float64, @inferred(logpdf(Normal(0.0, 0.0), NaN)))
    @test @inferred(logpdf(Normal(0.0f0, 0.0f0), 0.0))       === Inf
    @test @inferred(logpdf(Normal(0.0f0, 0.0f0), 0.0f0))     === Inf32
    @test isnan_type(Float64, @inferred(logpdf(Normal(0.0, 0.0), NaN)))
    @test isnan_type(Float64, @inferred(logpdf(Normal(NaN, 0.0), 0.0f0)))
    @test isnan_type(Float32, @inferred(logpdf(Normal(NaN32, 0.0f0), 0.0f0)))
    @test @inferred(logpdf(Normal(0 // 1, 0 // 1), 0 // 1))  === Inf
    @test isnan_type(Float64, @inferred(logpdf(Normal(0 // 1, 0 // 1), NaN)))
    @test @inferred(logpdf(Normal(0.0, 0.0), BigInt(1)))     == big(-Inf)
    @test @inferred(logpdf(Normal(0.0, 0.0), BigFloat(1)))   == big(-Inf)
    @test isnan_type(BigFloat, @inferred(logpdf(Normal(0.0, 0.0), BigFloat(NaN))))

    @test @inferred(cdf(Normal(0.0, 0.0), 0.0))           === 1.0
    @test @inferred(cdf(Normal(0.0, 0.0), -1.0))          === 0.0
    @test @inferred(cdf(Normal(0.0, 0.0), 0.0f0))         === 1.0
    @test isnan_type(Float64, @inferred(cdf(Normal(0.0, 0.0), NaN)))
    @test @inferred(cdf(Normal(0.0f0, 0.0f0), 0.0))       === 1.0
    @test @inferred(cdf(Normal(0.0f0, 0.0f0), 0.0f0))     === 1.0f0
    @test isnan_type(Float64, @inferred(cdf(Normal(0.0, 0.0), NaN)))
    @test isnan_type(Float64, @inferred(cdf(Normal(NaN, 0.0), 0.0f0)))
    @test isnan_type(Float32, @inferred(cdf(Normal(NaN32, 0.0f0), 0.0f0)))
    @test @inferred(cdf(Normal(0 // 1, 0 // 1), 0 // 1))  === 1.0
    @test isnan_type(Float64, @inferred(cdf(Normal(0 // 1, 0 // 1), NaN)))
    @test @inferred(cdf(Normal(0.0, 0.0), BigInt(1)))     == big(1.0)
    @test @inferred(cdf(Normal(0.0, 0.0), BigFloat(1)))   == big(1.0)
    @test isnan_type(BigFloat, @inferred(cdf(Normal(0.0, 0.0), BigFloat(NaN))))

    @test @inferred(logcdf(Normal(0.0, 0.0), 0.0))           === -0.0
    @test @inferred(logcdf(Normal(0.0, 0.0), -1.0))          === -Inf
    @test @inferred(logcdf(Normal(0.0, 0.0), 0.0f0))         === -0.0
    @test isnan_type(Float64, @inferred(logcdf(Normal(0.0, 0.0), NaN)))
    @test @inferred(logcdf(Normal(0.0f0, 0.0f0), 0.0))       === -0.0
    @test @inferred(logcdf(Normal(0.0f0, 0.0f0), 0.0f0))     === -0.0f0
    @test isnan_type(Float64, @inferred(logcdf(Normal(0.0, 0.0), NaN)))
    @test isnan_type(Float64, @inferred(logcdf(Normal(NaN, 0.0), 0.0f0)))
    @test isnan_type(Float32, @inferred(logcdf(Normal(NaN32, 0.0f0), 0.0f0)))
    @test @inferred(logcdf(Normal(0 // 1, 0 // 1), 0 // 1))  === -0.0
    @test isnan_type(Float64, @inferred(logcdf(Normal(0 // 1, 0 // 1), NaN)))
    @test @inferred(logcdf(Normal(0.0, 0.0), BigInt(1)))     == big(0.0)
    @test @inferred(logcdf(Normal(0.0, 0.0), BigFloat(1)))   == big(0.0)
    @test isnan_type(BigFloat, @inferred(logcdf(Normal(0.0, 0.0), BigFloat(NaN))))

    @test @inferred(ccdf(Normal(0.0, 0.0), 0.0))           === 0.0
    @test @inferred(ccdf(Normal(0.0, 0.0), -1.0))          === 1.0
    @test @inferred(ccdf(Normal(0.0, 0.0), 0.0f0))         === 0.0
    @test isnan_type(Float64, @inferred(ccdf(Normal(0.0, 0.0), NaN)))
    @test @inferred(ccdf(Normal(0.0f0, 0.0f0), 0.0))       === 0.0
    @test @inferred(ccdf(Normal(0.0f0, 0.0f0), 0.0f0))     === 0.0f0
    @test isnan_type(Float64, @inferred(ccdf(Normal(0.0, 0.0), NaN)))
    @test isnan_type(Float64, @inferred(ccdf(Normal(NaN, 0.0), 0.0f0)))
    @test isnan_type(Float32, @inferred(ccdf(Normal(NaN32, 0.0f0), 0.0f0)))
    @test @inferred(ccdf(Normal(0 // 1, 0 // 1), 0 // 1))  === 0.0
    @test isnan_type(Float64, @inferred(ccdf(Normal(0 // 1, 0 // 1), NaN)))
    @test @inferred(ccdf(Normal(0.0, 0.0), BigInt(1)))     == big(0.0)
    @test @inferred(ccdf(Normal(0.0, 0.0), BigFloat(1)))   == big(0.0)
    @test isnan_type(BigFloat, @inferred(ccdf(Normal(0.0, 0.0), BigFloat(NaN))))

    @test @inferred(logccdf(Normal(0.0, 0.0), 0.0))           === -Inf
    @test @inferred(logccdf(Normal(0.0, 0.0), -1.0))          === -0.0
    @test @inferred(logccdf(Normal(0.0, 0.0), 0.0f0))         === -Inf
    @test isnan_type(Float64, @inferred(logccdf(Normal(0.0, 0.0), NaN)))
    @test @inferred(logccdf(Normal(0.0f0, 0.0f0), 0.0))       === -Inf
    @test @inferred(logccdf(Normal(0.0f0, 0.0f0), 0.0f0))     === -Inf32
    @test isnan_type(Float64, @inferred(logccdf(Normal(0.0, 0.0), NaN)))
    @test isnan_type(Float64, @inferred(logccdf(Normal(NaN, 0.0), 0.0f0)))
    @test isnan_type(Float32, @inferred(logccdf(Normal(NaN32, 0.0f0), 0.0f0)))
    @test @inferred(logccdf(Normal(0 // 1, 0 // 1), 0 // 1))  === -Inf
    @test isnan_type(Float64, @inferred(logccdf(Normal(0 // 1, 0 // 1), NaN)))
    @test @inferred(logccdf(Normal(0.0, 0.0), BigInt(1)))     == big(-Inf)
    @test @inferred(logccdf(Normal(0.0, 0.0), BigFloat(1)))   == big(-Inf)
    @test isnan_type(BigFloat, @inferred(logccdf(Normal(0.0, 0.0), BigFloat(NaN))))

    @test @inferred(quantile(Normal(1.0, 0.0), 0.0f0))     === -Inf
    @test @inferred(quantile(Normal(1.0, 0.0f0), 1.0))     ===  Inf
    @test @inferred(quantile(Normal(1.0f0, 0.0), 0.5))     ===  1.0
    @test isnan_type(Float64, @inferred(quantile(Normal(1.0f0, 0.0), NaN)))
    @test @inferred(quantile(Normal(1.0f0, 0.0f0), 0.0f0)) === -Inf32
    @test @inferred(quantile(Normal(1.0f0, 0.0f0), 1.0f0)) ===  Inf32
    @test @inferred(quantile(Normal(1.0f0, 0.0f0), 0.5f0)) ===  1.0f0
    @test isnan_type(Float32, @inferred(quantile(Normal(1.0f0, 0.0f0), NaN32)))
    @test @inferred(quantile(Normal(1//1, 0//1), 1//2))    ===  1.0

    @test @inferred(cquantile(Normal(1.0, 0.0), 0.0f0))     ===  Inf
    @test @inferred(cquantile(Normal(1.0, 0.0f0), 1.0))     === -Inf
    @test @inferred(cquantile(Normal(1.0f0, 0.0), 0.5))     ===  1.0
    @test isnan_type(Float64, @inferred(cquantile(Normal(1.0f0, 0.0), NaN)))
    @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 0.0f0)) ===  Inf32
    @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 1.0f0)) === -Inf32
    @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 0.5f0)) ===  1.0f0
    @test isnan_type(Float32, @inferred(cquantile(Normal(1.0f0, 0.0f0), NaN32)))
    @test @inferred(cquantile(Normal(1//1, 0//1), 1//2))    ===  1.0
end

@testset "Normal: Sampling with integer-valued parameters" begin
    d = Normal{Int}(0, 1)
    @test rand(d) isa Float64
    @test rand(d, 10) isa Vector{Float64}
    @test rand(d, (3, 2)) isa Matrix{Float64}
end

@testset "Normal fit to stats" begin
    d = Normal(3,2)
    qp = @qp(quantile(d,0.95),0.95)
    @testset "fit moments" begin
        dfit = fit(Normal, moments(d))
        @test dfit ≈ d
        @test_throws ErrorException d = fit(Normal, Moments(3.0))
    end;
    @testset "fit two quantiles" begin
        qpl = @qp_m(3)
        qpu = @qp_u(5)
        dn = fit(Normal, qpl, qpu)
        @test quantile.(dn, [qpl.p, qpu.p]) ≈ [qpl.q, qpu.q]
        dn = fit(Normal, qpu, qpl) # check resorting arguments
        @test quantile.(dn, [qpl.p, qpu.p]) ≈ [qpl.q, qpu.q]
    end;
    @testset "fit to quantilepoint and mean" begin
        dfit = fit(Normal, mean(d), qp, Val(:mean))
        @test dfit ≈ d
    end;
    @testset "fit to quantilepoint and mode" begin
        dfit = fit(Normal, mode(d), qp, Val(:mode))
        @test dfit ≈ d
    end;
    @testset "fit to quantilepoint and median" begin
        dfit = fit(Normal, median(d), qp, Val(:median))
        @test dfit ≈ d
    end;
end;