text:
export dataset,
       findVar,
       isArray,
       isFactor,
       isMatrix,
       lang1,
       lang2,
       lang3,
       lang4,
       lang5,
       lang6,
       reval,
       rparse,
       rprint,
       sexp

@doc "evaluate an R symbol or language object (i.e. a function call) in an R try/catch block"->
function reval(expr::SEXP, env::SEXP{ENVSXP})
    err = Array(Cint,1)
    val = ccall((:R_tryEval,libR),Ptr{Void},(Ptr{Void},Ptr{Void},Ptr{Cint}),expr,env,err)
    bool(err[1]) && error("Error occurred in R_tryEval")
    sexp(val)
end

@doc "expression objects (the result of rparse) have a special reval method"->
function reval(expr::SEXP{EXPRSXP}, env::SEXP{ENVSXP}) # evaluate result of R_ParseVector
    local val           # the value of the last expression is returned
    for e in expr
        val = reval(e,env)
    end
    val
end
reval(s::SEXP) = reval(s,globalEnv)
reval(sym::Symbol) = reval(sexp(sym))

@doc "return the first element of an SEXP as an Complex128 value" ->
asComplex(s::SEXP) = ccall((:Rf_asComplex,libR),Complex128,(Ptr{Void},),s)

@doc "return the first element of an SEXP as an Cint (i.e. Int32)" ->
asInteger(s::SEXP) = ccall((:Rf_asInteger,libR),Cint,(Ptr{Void},),s)

@doc "return the first element of an SEXP as a Bool" ->
asLogical(s::SEXP) = ccall((:Rf_asLogical,libR),Bool,(Ptr{Void},),s)

@doc "return the first element of an SEXP as a Cdouble (i.e. Float64)" ->
asReal(s::SEXP) = ccall((:Rf_asReal,libR),Cdouble,(Ptr{Void},),s)

@doc "find object with name sym in environment env"->
findVar(sym::SEXP,env::SEXP{ENVSXP}=globalEnv) =
    sexp(ccall((:Rf_findVar,libR),Ptr{Void},(Ptr{Void},Ptr{Void}),sym,env))
findVar(nm::ASCIIString,env::SEXP{ENVSXP}) = findVar(sexp(symbol(nm)),env)
findVar(nm::ASCIIString) = findVar(sexp(symbol(nm)),globalEnv)

## predicates applied to an SEXP (many of these are unneeded for templated SEXP)
for sym in (:isArray,:isComplex,:isEnvironment,:isExpression,:isFactor,
            :isFrame,:isFree,:isFunction,:isInteger,:isLanguage,:isList,
            :isLogical,:isSymbol,:isMatrix,:isNewList,:isNull,:isNumeric,
            :isNumber,:isObject,:isOrdered,:isPairList,:isPrimitive,
            :isReal,:isS4,:isString,:isTs,:isUnordered,:isUnsorted,
            :isUserBinop,:isValidString,:isValidStringF,:isVector,
            :isVectorAtomic,:isVectorizable,:isVectorList)
    @eval $sym(s::SEXP) = ccall(($(string("Rf_",sym)),libR),Bool,(Ptr{Void},),s)
end

@doc "Create a 0-argument function call from a symbol"->
lang1(s::SEXP) = sexp(ccall((:Rf_lang1,libR),Ptr{Void},(Ptr{Void},),s))

@doc "Create a 1-argument function call from a symbol and the argument"->
lang2(s1::SEXP,s2::SEXP) =
    sexp(ccall((:Rf_lang2,libR),Ptr{Void},(Ptr{Void},Ptr{Void}),s1,s2))

@doc "Create a 2-argument function call from a symbol and the arguments"->
lang3(s1::SEXP,s2::SEXP,s3::SEXP) =
    sexp(ccall((:Rf_lang3,libR),Ptr{Void},
                 (Ptr{Void},Ptr{Void},Ptr{Void}),s1,s2,s3))

@doc "Create a 3-argument function call from a symbol and the arguments"->
lang4(s1::SEXP,s2::SEXP,s3::SEXP,s4::SEXP) =
    sexp(ccall((:Rf_lang3,libR),Ptr{Void},
                 (Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),s1,s2,s3,s4))

@doc "Create a 4-argument function call from a symbol and the arguments"->
lang5(s1::SEXP,s2::SEXP,s3::SEXP,s4::SEXP,s5::SEXP) =
    sexp(ccall((:Rf_lang3,libR),Ptr{Void},
                 (Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),s1,s2,s3,s4,s5))

@doc "Create a 5-argument function call from a symbol and the arguments"->
lang6(s1::SEXP,s2::SEXP,s3::SEXP,s4::SEXP,s5::SEXP,s6::SEXP) =
    sexp(ccall((:Rf_lang3,libR),Ptr{Void},
                 (Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void},Ptr{Void}),
                 s1,s2,s3,s4,s5,s6))

@doc "Protect an SEXP from garbage collection"->
protect(s::SEXP) = sexp(ccall((:Rf_protect,libR),Ptr{Void},(Ptr{Void},),s))

@doc "Parse a string as an R expression"->
function rparse(st::ASCIIString)
    ParseStatus = Array(Cint,1)
    val = ccall((:R_ParseVector,libR),Ptr{Void},
                (Ptr{Void},Cint,Ptr{Cint},Ptr{Void}),
                sexp(st),length(st),ParseStatus,nilValue)
    ParseStatus[1] == 1 || error("R_ParseVector set ParseStatus to $(ParseStatus[1])")
    sexp(val)
end

@doc "print the value of an SEXP using R's printing mechanism"->
rprint(s::SEXP) = ccall((:Rf_PrintValue,libR),Void,(Ptr{Void},),s)

@doc "Pop k elements off the protection stack"->
unprotect(k::Integer) = ccall((:Rf_unprotect,libR),Void,(Cint,),k)

@doc "unprotect an SEXP"->
unprotect(s::SEXP) = ccall((:Rf_unprotect_ptr,libR),Void,(Ptr{Void},),s)

@doc "release an SEXP"->
ReleaseObject(s::SEXP) = ccall((:R_ReleaseObject,libR),Void,(Ptr{Void},),s)

@doc "preserve an SEXP"->
PreserveObject(s::SEXP) = ccall((:R_PreserveObject,libR),Void,(Ptr{Void},),s)


text:
@testset "Quadrature testing" begin
    ref_tet_vol(dim) = 1 / factorial(dim)
    ref_square_vol(dim) = 2^dim

    function integrate(qr::QuadratureRule, f::Function)
        I = 0.0
        for (w, x) in zip(getweights(qr), getpoints(qr))
            I += w * f(x)
        end
        return I
    end

    # Cube
    for dim = (1,2,3)
        for order in (1,2,3,4)
            f = (x, p) -> sum([x[i]^p for i in 1:length(x)])
            # Legendre
            qr = QuadratureRule{dim, RefCube}(:legendre, order)
            @test integrate(qr, (x) -> f(x, 2*order-1)) < 1e-14
            @test sum(qr.weights) ≈ ref_square_vol(dim)
            @test sum(getweights(qr)) ≈ ref_square_vol(dim)
            # Lobatto
            if order > 1
                qr = QuadratureRule{dim, RefCube}(:lobatto, order)
                @test integrate(qr, (x) -> f(x, 2*order-1)) < 1e-14
                @test sum(qr.weights) ≈ ref_square_vol(dim)
                @test sum(getweights(qr)) ≈ ref_square_vol(dim)
            end
        end
    end
    @test_throws ArgumentError QuadratureRule{1, RefCube}(:einstein, 2)

    # Tetrahedron
    g = (x) -> sqrt(sum(x))
    dim = 2
    for order in (1,2,3,4,5,6,7,8)
        qr = QuadratureRule{dim, RefTetrahedron}(:legendre, order)
        # http://www.wolframalpha.com/input/?i=integrate+sqrt(x%2By)+from+x+%3D+0+to+1,+y+%3D+0+to+1-x
        @test integrate(qr, g) - 0.4 < 0.01
        @test sum(qr.weights) ≈ ref_tet_vol(dim)
    end
    @test_throws ArgumentError QuadratureRule{dim, RefTetrahedron}(:einstein, 2)
    @test_throws ArgumentError QuadratureRule{dim, RefTetrahedron}(0)

    dim = 3
    for order in (1, 2, 3)
        qr = QuadratureRule{dim, RefTetrahedron}(:legendre, order)
        # Table 1:
        # http://www.m-hikari.com/ijma/ijma-2011/ijma-1-4-2011/venkateshIJMA1-4-2011.pdf
        @test integrate(qr, g) - 0.14 < 0.01
        @test sum(qr.weights) ≈ ref_tet_vol(dim)
    end
    @test_throws ArgumentError QuadratureRule{dim, RefTetrahedron}(:einstein, 2)
    @test_throws ArgumentError QuadratureRule{dim, RefTetrahedron}(0)

end

ext: jl


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

# math functions
libm_math_functions = Set([:sin, :cos, :tan, :asin, :acos, :acosh, :atanh, :log, :log2, :log10, :lgamma, :log1p,:asinh,:atan,:cbrt,:cosh,:erf,:exp,:expm1,:sinh,:sqrt,:tanh, :isnan])
#using Debug

function pattern_match_call_math(fun::Symbol, input::AbstractString, typ::Type, linfo)
    s = ""
    isDouble = typ == Float64
    isFloat = typ == Float32
    isComplex = typ <: Complex
    isInt = typ <: Integer
    if in(fun,libm_math_functions) && (isFloat || isDouble || isComplex)
        @dprintln(3,"FOUND ", fun)
        s = string(fun)*"("*input*");"
    end

    # abs() needs special handling since fabs() in math.h should be called for floats
    if is(fun,:abs) && (isFloat || isDouble || isComplex || isInt)
      @dprintln(3,"FOUND ", fun)
      fname = (isInt || isComplex) ? "abs" : (isFloat ? "fabsf" : "fabs")
      s = fname*"("*input*");"
    end
    return s
end

function pattern_match_call_math(fun::Symbol, input::RHSVar, linfo)
  pattern_match_call_math(fun, from_expr(input, linfo), getType(input, linfo), linfo)
end

function pattern_match_call_math(fun::GlobalRef, input, linfo)
    fun = Base.resolve(fun)
    if fun.mod == Base
        pattern_match_call_math(fun.name, input,linfo)
    else
        return ""
    end
end

function pattern_match_call_math(fun::ANY, input::ANY, linfo)
    return ""
end

function pattern_match_call_throw(fun::GlobalRef, input, linfo)
    s = ""
    if fun.name==:throw || fun.name==:error
        s = "(throw(\"Julia throw() or error() called.\"), 0)"
    end
    return s
end

function pattern_match_call_throw(fun::Symbol, input, linfo)
    s = ""
    if fun==:throw || fun==:error
        s = "(throw(\"Julia throw() or error() called.\"), 0)"
    end
    return s
end

function pattern_match_call_throw(fun::ANY, input::ANY, linfo)
    return ""
end

function pattern_match_call_powersq(fun, x::Number, y::Integer, linfo)
    s = ""
    if isBaseFunc(fun, :power_by_squaring)
        s = "cgen_pown("*from_expr(x,linfo)*","*from_expr(y,linfo)*")"
    end
    return s
end

function pattern_match_call_powersq(fun::ANY, x::ANY, y::ANY,linfo)
    return ""
end

function pattern_match_call_rand(linfo, fun, args...)
    @dprintln(3,"pattern_match_call_rand ", fun)
    res = ""
    if isBaseFunc(fun, :rand)
        if USE_OMP==1
            res = "cgen_distribution(cgen_rand_generator[omp_get_thread_num()]);\n"
        else
            res = "cgen_distribution(cgen_rand_generator);\n"
        end
    end
    @dprintln(3,"pattern_match_call_rand res = ", res)
    return res
end

function pattern_match_call_randn(linfo, fun, args...)
    @dprintln(3,"pattern_match_call_randn ", fun)
    res = ""
    if isBaseFunc(fun, :randn)
        if USE_OMP==1
            res = "cgen_n_distribution(cgen_rand_generator[omp_get_thread_num()]);\n"
        else
            res = "cgen_n_distribution(cgen_rand_generator);\n"
        end
    end
    @dprintln(3,"pattern_match_call_randn res = ", res)
    return res
end

function pattern_match_call_reshape(fun, inp::Any, shape::RHSVar, linfo)
    res = ""
    if isBaseFunc(fun, :reshape)
        typ = getSymType(shape, linfo)
        if istupletyp(typ)
            dim = length(typ.parameters)
            sh = from_expr(shape,linfo)
            shapes = mapfoldl(i->sh*".f"*string(i-1), (a,b) -> a*","*b, 1:dim)
            res = from_expr(inp,linfo) * ".reshape(" * shapes * ");\n"
        else
            error("call to reshape expects a tuple, but got ", typ)
        end
    end
    return res
end

function pattern_match_call_reshape(fun::ANY, inp::ANY, shape::ANY,linfo)
    return ""
end

function getSymType(a, linfo)
    return lstate.symboltable[lookupVariableName(a, linfo)]
end

function pattern_match_call_gemm(fun::GlobalRef, C::RHSVar, tA::Char, tB::Char, A::RHSVar, B::RHSVar,linfo)
    if fun.mod!=Base.LinAlg || fun.name!=:gemm_wrapper!
        return ""
    end
    cblas_fun = ""
    typ = getSymType(A, linfo)
    if getSymType(B, linfo) != typ || getSymType(C, linfo) != typ
        return ""
    end
    if typ==Array{Float32,2}
        cblas_fun = "cblas_sgemm"
    elseif typ==Array{Float64,2}
        cblas_fun = "cblas_dgemm"
    else
        return ""
    end
    s = "$(from_expr(C,linfo)); "
    # GEMM wants dimensions after possible transpose
    m = (tA == 'N') ? from_arraysize(A,1,linfo) : from_arraysize(A,2,linfo)
    k = (tA == 'N') ? from_arraysize(A,2,linfo) : from_arraysize(A,1,linfo)
    n = (tB == 'N') ? from_arraysize(B,2,linfo) : from_arraysize(B,1,linfo)

    lda = from_arraysize(A,1,linfo)
    ldb = from_arraysize(B,1,linfo)
    ldc = m

    CblasNoTrans = 111
    CblasTrans = 112
    _tA = tA == 'N' ? CblasNoTrans : CblasTrans
    _tB = tB == 'N' ? CblasNoTrans : CblasTrans
    CblasColMajor = 102

    if mkl_lib!="" || openblas_lib!="" || sys_blas==1
        s *= "$(cblas_fun)((CBLAS_ORDER)$(CblasColMajor),(CBLAS_TRANSPOSE)$(_tA),(CBLAS_TRANSPOSE)$(_tB),$m,$n,$k,1.0,
        $(from_expr(A,linfo)).data, $lda, $(from_expr(B,linfo)).data, $ldb, 0.0, $(from_expr(C,linfo)).data, $ldc)"
    else
        println("WARNING: MKL and OpenBLAS not found. Matrix multiplication might be slow.
        Please install MKL or OpenBLAS and rebuild ParallelAccelerator for better performance.")
        s *= "cgen_$(cblas_fun)($(from_expr(tA!='N',linfo)), $(from_expr(tB!='N',linfo)), $m,$n,$k, $(from_expr(A,linfo)).data, $lda, $(from_expr(B,linfo)).data, $ldb, $(from_expr(C,linfo)).data, $ldc)"
    end

    return s
end

function pattern_match_call_gemm(fun::ANY, C::ANY, tA::ANY, tB::ANY, A::ANY, B::ANY,linfo)
    return ""
end

function pattern_match_call_gemv(fun::GlobalRef, y::RHSVar, tA::Char, A::RHSVar, x::RHSVar,linfo)
    if !((fun.mod==Base.LinAlg || fun.mod==Base.LinAlg.BLAS) && fun.name==:gemv!)
        return ""
    end
    cblas_fun = ""
    typ = eltype(getSymType(A, linfo))

    if typ==Float32
        cblas_fun = "cblas_sgemv"
    elseif typ==Float64
        cblas_fun = "cblas_dgemv"
    else
        return ""
    end

    s = "$(from_expr(y,linfo)); "

    m = from_arraysize(A,1,linfo)
    n = from_arraysize(A,2,linfo)

    lda = from_arraysize(A,1,linfo)

    CblasNoTrans = 111
    CblasTrans = 112
    _tA = tA == 'N' ? CblasNoTrans : CblasTrans
    CblasColMajor = 102

    if mkl_lib!="" || openblas_lib!="" || sys_blas==1
        s *= "$(cblas_fun)((CBLAS_ORDER)$(CblasColMajor),(CBLAS_TRANSPOSE)$(_tA),$m,$n, 1.0,
        $(from_expr(A,linfo)).data, $lda, $(from_expr(x,linfo)).data, 1, 0.0, $(from_expr(y,linfo)).data, 1)"
    else
        println("WARNING: MKL and OpenBLAS not found. Matrix-vector multiplication might be slow.
        Please install MKL or OpenBLAS and rebuild ParallelAccelerator for better performance.")
        s *= "cgen_$(cblas_fun)($(from_expr(tA!='N',linfo)), $m,$n, $(from_expr(A,linfo)).data, $lda, $(from_expr(y,linfo)).data, $(from_expr(x,linfo)).data)"
    end

    return s
end

function pattern_match_call_gemv(fun::ANY, C::ANY, tA::ANY, A::ANY, B::ANY,linfo)
    return ""
end

function pattern_match_call_chol(fun::GlobalRef, A::RHSVar, vUL::Type, linfo)
    if fun.mod!=Base.LinAlg || fun.name!=:chol!
        return ""
    end

    cblas_fun = ""
    typ = eltype(getSymType(A, linfo))

    if typ==Float32
        lapack_fun = "LAPACKE_spotrf"
    elseif typ==Float64
        lapack_fun = "LAPACKE_dpotrf"
    else
        return ""
    end

    s = ".data=$(from_expr(A,linfo)); "

    n = from_arraysize(A,1,linfo)

    lda = from_arraysize(A,1,linfo)

    LAPACK_COL_MAJOR = 102
    uplo = vUL==Val{:U} ? 'U' : 'L'

    if mkl_lib!="" || openblas_lib!="" || sys_blas==1
        s *= "$(lapack_fun)($(LAPACK_COL_MAJOR), '$uplo', $n, $(from_expr(A,linfo)).data, $lda)"
    else
        println("WARNING: MKL and OpenBLAS not found. Matrix multiplication might be slow.
        Please install MKL or OpenBLAS and rebuild ParallelAccelerator for better performance.")
        error("MKL LAPACK required for cholesky (TODO: support other lapack libraries and include a sequential implementation)")
        #s *= "cgen_$(cblas_fun)($(from_expr(tA!='N',linfo)), $m,$n, $(from_expr(A,linfo)).data, $lda, $(from_expr(y,linfo)).data, $(from_expr(x,linfo)).data)"
    end

    return s
end

function pattern_match_call_chol(fun::ANY, C::ANY, tA::ANY, linfo)
    return ""
end

function pattern_match_assignment_chol(lhs::LHSVar, rhs::Expr, linfo)
    call = ""
    if isCall(rhs) || isInvoke(rhs)
        fun = getCallFunction(rhs)
        args = getCallArguments(rhs)
        if length(args) == 2
            call *= pattern_match_call_chol(fun,args[1],args[2],linfo)
        end
    end
    if call!=""
        return from_expr(lhs,linfo)*call
    end
    return ""
end

function pattern_match_assignment_chol(lhs::ANY, rhs::ANY, linfo)
    return ""
end

function pattern_match_assignment_transpose(lhs::LHSVar, rhs::Expr, linfo)
    @dprintln(3, "pattern_match_assignment_transpose ", lhs, " ", rhs)
    call = ""
    if isCall(rhs) || isInvoke(rhs)
        fun = getCallFunction(rhs)
        args = getCallArguments(rhs)
        res = pattern_match_call_transpose(linfo, fun, lhs, args...)
        return res
    end
    return ""
end

function pattern_match_assignment_transpose(lhs::ANY, rhs::ANY, linfo)
    return ""
end

function pattern_match_call_trtrs(fun::GlobalRef, uplo::Char, trans::Char, diag::Char, A::RHSVar,B::RHSVar,  linfo)
    if fun.mod!=Base.LinAlg.LAPACK || fun.name!=:trtrs!
        return ""
    end

    cblas_fun = ""
    typ = eltype(getSymType(A, linfo))

    if typ==Float32
        lapack_fun = "LAPACKE_strtrs"
    elseif typ==Float64
        lapack_fun = "LAPACKE_dtrtrs"
    else
        return ""
    end

    s = "$(from_expr(B,linfo)); "

    n = from_arraysize(A,1,linfo)
    nrhs = from_arraysize(B,2,linfo)

    lda = from_arraysize(A,1,linfo)
    ldb = n

    LAPACK_COL_MAJOR = 102

    if mkl_lib!="" || openblas_lib!="" || sys_blas==1
        s *= "$(lapack_fun)($(LAPACK_COL_MAJOR), '$uplo', '$trans', '$diag', $n, $nrhs, $(from_expr(A,linfo)).data, $lda, $(from_expr(B,linfo)).data, $ldb)"
    else
        println("WARNING: MKL and OpenBLAS not found. Matrix multiplication might be slow.
        Please install MKL or OpenBLAS and rebuild ParallelAccelerator for better performance.")
        error("MKL LAPACK required for triangular solution (TODO: support other lapack libraries and include a sequential implementation)")
        #s *= "cgen_$(cblas_fun)($(from_expr(tA!='N',linfo)), $m,$n, $(from_expr(A,linfo)).data, $lda, $(from_expr(y,linfo)).data, $(from_expr(x,linfo)).data)"
    end

    return s
end

function pattern_match_call_trtrs(fun::ANY, C::ANY, tA::ANY,A::ANY, t::ANY, d::ANY, linfo)
    return ""
end

function pattern_match_call_copy!(linfo, fun, A, i, B, j, n)
    if isBaseFunc(fun, :copy!)
        "j2c_array_copyto(" *
          from_expr(A, linfo) * ", " *
          from_expr(i, linfo) * ", " *
          from_expr(B, linfo) * ", " *
          from_expr(j, linfo) * ", " *
          from_expr(n, linfo) * ")"
    else
        ""
    end
end

function pattern_match_call_transpose(linfo, fun::GlobalRef, fun1::GlobalRef, B::RHSVar, A::RHSVar)
    pattern_match_call_transpose(linfo, fun, B, A)
end

function pattern_match_call_transpose(linfo, fun::GlobalRef, B::RHSVar, A::RHSVar)
    dprintln(3, "pattern_match_call_transpose, ", (fun, B, A), " mkl_lib=", mkl_lib, " openblas=",openblas_lib, " sys_blas=",sys_blas)
    if !(fun.mod==Base && fun.name==:transpose! || fun.name ==:transpose_f! || fun.name == :transpose)
        return ""
    end
    blas_fun = ""
    typ = eltype(getSymType(A, linfo))
    ctyp = toCtype(typ)

    if typ==Float32
        blas_fun = "somatcopy"
    elseif typ==Float64
        blas_fun = "domatcopy"
    else
        blas_fun = ""
    end

    #s = "$(from_expr(B,linfo)); "
    s = ""

    m = from_arraysize(A,1,linfo)
    n = from_arraysize(A,2,linfo)

    lda = from_arraysize(A,1,linfo)
    ldb = from_arraysize(A,2,linfo)

    if mkl_lib!="" && blas_fun!=""
        s *= "mkl_$(blas_fun)('C','T',$m,$n, 1.0,
             $(from_expr(A,linfo)).data, $lda, $(from_expr(B,linfo)).data, $ldb)"
    elseif (openblas_lib!="" || sys_blas==1) && blas_fun!=""
        s *= "cblas_$(blas_fun)(CblasColMajor,CblasTrans,$m,$n, 1.0,
             $(from_expr(A,linfo)).data, $lda, $(from_expr(B,linfo)).data, $ldb)"
    else
        #println("""WARNING: MKL and OpenBLAS not found. Matrix-vector multiplication might be slow.
        #Please install MKL or OpenBLAS and rebuild ParallelAccelerator for better performance.""")
        #s *= "cgen_$(blas_fun)($m,$n,
        #     $(from_expr(A,linfo)).data, $lda, $(from_expr(B,linfo)).data, $ldb)"
        s *= "for(int i=0; i<$m; i++) {\n"
        s *= "    for(int j=0; j<$n; j++) {\n"
        s *= "     $(from_expr(B,linfo)).data[j+i*$ldb] = $(from_expr(A,linfo)).data[i+j*$lda];\n"
        s *= "    }\n"
        s *= "}\n"
    end
    s = (fun.name == :transpose ? (from_expr(B,linfo) * " = j2c_array<$ctyp>::new_j2c_array_2d(NULL, $n, $m);\n"): "") * s
    s = from_expr(B,linfo)*"; "*s
    return s
end

function pattern_match_call_transpose(args...)
    return ""
end

function pattern_match_call_linalgtypeof(fun::GlobalRef, C::ANY,linfo)
    if fun.mod==Base.LinAlg && fun.name==:typeof
        return " "
    end
    return ""
end

function pattern_match_call_linalgtypeof(fun::ANY, C::ANY,linfo)
    return ""
end

function pattern_match_call_vecnorm(fun::GlobalRef, y::RHSVar, p::Int,linfo)
    if fun.mod!=Base.LinAlg || fun.name!=:vecnorm
        return ""
    end
    cblas_fun = ""
    typ = eltype(getSymType(y, linfo))

    if typ==Float32
        cblas_fun = "cblas_s"
    elseif typ==Float64
        cblas_fun = "cblas_d"
    else
        return ""
    end

    if p==1
        cblas_fun *= "asum"
    elseif p==2
        cblas_fun *= "nrm2"
    else
        println("norm ",p)
        error("vector norm not support")
    end

    s = ""

    n = from_arraysize(y,1,linfo)

    if mkl_lib!="" || openblas_lib!="" || sys_blas==1
        s *= "$(cblas_fun)( $n, $(from_expr(y,linfo)).data, 1)"
    else
        #println("WARNING: MKL and OpenBLAS not found. Matrix-vector multiplication might be slow.
        #Please install MKL or OpenBLAS and rebuild ParallelAccelerator for better performance.")
        s *= "cgen_$(cblas_fun)( $n, $(from_expr(y,linfo)).data)"
    end

    return s
end

function pattern_match_call_vecnorm(fun::ANY, C::ANY, tA::ANY,linfo)
    return ""
end

function pattern_match_call(ast::Array{Any, 1},linfo)
    @dprintln(3,"pattern matching ",ast)
    s = ""

    if(length(ast)==2)
        s = pattern_match_call_throw(ast[1],ast[2],linfo)
        s *= pattern_match_call_math(ast[1],ast[2],linfo)
        s *= pattern_match_call_linalgtypeof(ast[1],ast[2],linfo)
    end

    if s=="" && (length(ast)==3) # randn! call has 3 args
        #sa*= pattern_match_call_powersq(ast[1],ast[2], ast[3])
        s *= pattern_match_call_reshape(ast[1],ast[2],ast[3],linfo)
        s *= pattern_match_call_transpose(linfo, ast...)
        s *= pattern_match_call_vecnorm(ast[1],ast[2],ast[3],linfo)
    end
    if s=="" && (length(ast)>=1) # rand can have 1 or more arg
        s *= pattern_match_call_transpose(linfo, ast...)
        s *= pattern_match_call_randn(linfo, ast...)
        s *= pattern_match_call_rand(linfo, ast...)
    end
    # gemv calls have 5 args
    if s=="" && (length(ast)==5)
        s *= pattern_match_call_gemv(ast[1],ast[2],ast[3],ast[4],ast[5],linfo)
    end
    # gemm calls have 6 args
    if s=="" && (length(ast)==6)
        s *= pattern_match_call_copy!(linfo, ast...)
        s *= pattern_match_call_gemm(ast[1],ast[2],ast[3],ast[4],ast[5],ast[6],linfo)
        s *= pattern_match_call_trtrs(ast[1],ast[2],ast[3],ast[4],ast[5],ast[6],linfo)
    end
    return s
end

function from_assignment_match_hvcat(lhs, rhs::Expr, linfo)
    s = ""
    # if this is a hvcat call, the array should be allocated and initialized
    if (isCall(rhs) || isInvoke(rhs)) && (isBaseFunc(getCallFunction(rhs),:typed_hvcat) || checkGlobalRefName(getCallFunction(rhs),:hvcat))
        @dprintln(3,"Found hvcat assignment: ", lhs," ", rhs)

        is_typed::Bool = isBaseFunc(getCallFunction(rhs),:typed_hvcat)

        rows = Int64[]
        values = Any[]
        typ = "double"
        args = getCallArguments(rhs)

        if is_typed
            atyp = args[1]
            if isa(atyp, GlobalRef)
                atyp = eval(args[1].name)
            end
            @assert isa(atyp, DataType) ("hvcat expects the first argument to be a type, but got " * args[1])
            typ = toCtype(atyp)
            rows = lstate.tupleTable[args[2]]
            values = args[3:end]
        else

            rows = lstate.tupleTable[args[1]]
            values = args[2:end]
            atyp, arr_dims = parseArrayType(getSymType(lhs, linfo))
            typ = toCtype(atyp)
        end

        nr = length(rows)
        nc = rows[1] # all rows should have the same size
        s *= from_expr(lhs,linfo) * " = j2c_array<$typ>::new_j2c_array_2d(NULL, $nr, $nc);\n"
        s *= mapfoldl((i) -> from_setindex([lhs,values[i],convert(Int64,ceil(i/nc)),(i-1)%nc+1],linfo)*";", (a, b) -> "$a $b", 1:length(values))
    end
    return s
end

function from_assignment_match_hvcat(lhs, rhs::ANY, linfo)
    return ""
end

function from_assignment_match_cat_t(lhs, rhs::Expr, linfo)
    s = ""
    if (isCall(rhs) || isInvoke(rhs)) && isBaseFunc(getCallFunction(rhs), :cat_t)
        args = getCallArguments(rhs)
        dims = args[1]
        @assert dims==2 "CGen: only 2d cat_t() is supported now"
        size = length(args[3:end])
        typ = toCtype(eval(args[2].name))
        s *= from_expr(lhs,linfo) * " = j2c_array<$typ>::new_j2c_array_$(dims)d(NULL, 1,$size);\n"
        values = args[3:end]
        s *= mapfoldl((i) -> from_setindex([lhs,values[i],i],linfo)*";", (a, b) -> "$a $b", 1:length(values))
    end
    return s
end

function from_assignment_match_cat_t(lhs, rhs::ANY, linfo)
    return ""
end

function from_assignment_match_hcat(lhs, rhs::Expr, linfo)
    s = ""
    if (isCall(rhs) || isInvoke(rhs)) && isBaseFunc(getCallFunction(rhs), :hcat)
        args = getCallArguments(rhs)
        in_typ = getType(args[1], linfo)
        # hcat of single values like hcat(1,2,3) => 2D array [1 2 3]
        # used in quant example
        if !(in_typ<:Array)
            return single_value_hcat(lhs, args, linfo)
        end
        for a in args
            atyp = getType(a, linfo)
            @assert atyp<:Array && ndims(atyp)==1 "CGen only supports hcat of 1D arrays"
        end
        typ = eltype(getType(lhs, linfo))
        size = length(args)
        ctyp = toCtype(typ)
        len = from_arraysize(args[1],1,linfo)
        clhs = from_expr(lhs,linfo)
        s *= "$clhs = j2c_array<$ctyp>::new_j2c_array_2d(NULL, $len, $size);\n"
        s *= "for(int i=0; i<$len; i++) {\n"
        for j in 1:size
            arr = from_expr(args[j],linfo)
            s *= "$clhs.data[i*$size+$(j-1)] = $arr.data[i];\n"
        end
        s *= "}\n"
    end
    return s
end

function from_assignment_match_hcat(lhs, rhs::ANY, linfo)
    return ""
end

function single_value_hcat(lhs, args, linfo)
    for a in args
        atyp = getType(a, linfo)
        @assert !(atyp<:Array) "CGen invalid hcat input (single_value_hcat)"
    end
    typ = getType(args[1], linfo)
    size = length(args)
    ctyp = toCtype(typ)
    clhs = from_expr(lhs,linfo)
    s = "$clhs = j2c_array<$ctyp>::new_j2c_array_2d(NULL, 1, $size);\n"
    for j in 1:size
        var = from_expr(args[j],linfo)
        s *= "$clhs.data[$(j-1)] = $var;\n"
    end
    return s
end

function from_assignment_match_vcat(lhs, rhs::Expr, linfo)
    s = ""
    if (isCall(rhs) || isInvoke(rhs)) && isBaseFunc(getCallFunction(rhs), :vcat)
        args = getCallArguments(rhs)
        for a in args
            atyp = getType(a, linfo)
            @assert atyp<:Array && ndims(atyp)==1 "CGen only supports vcat of 1D arrays"
        end
        typ = eltype(getType(args[1], linfo))
        ctyp = toCtype(typ)
        clhs = from_expr(lhs,linfo)
        # get total size of array: size(a1)+size(a2)+...
        csize = "("* mapfoldl(a->from_arraysize(a,1,linfo),(a,b)->"$a+$b",args) *")"
        s *= "{\n"
        s *= "$clhs = j2c_array<$ctyp>::new_j2c_array_1d(NULL, $csize);\n"
        s *= "  int64_t __cgen_curr_ind = 0;\n"
        for arr in args
            carr = from_expr(arr, linfo)
            s *= "for(int64_t i=0; i<$(from_arraysize(arr,1,linfo)); i++){\n"
            s *= "  $clhs.data[__cgen_curr_ind++] = $carr.data[i];\n"
            s *= "}\n"
        end
        s *= "}\n"
    end
    return s
end

function from_assignment_match_vcat(lhs, rhs::ANY, linfo)
    return ""
end

function from_assignment_match_iostream(lhs, rhs::GlobalRef, linfo)
    s = ""
    ltype = getType(lhs, linfo)
    @dprintln(3, "from_assignment_match_iostream ltype = ", ltype)
    if (ltype == IOStream)
        if rhs.mod == Base && rhs.name == :STDOUT
            lhsO = from_expr(lhs, linfo)
            s *= lhsO * ".handle = (void**)&(std::cout);"
        end
    end
    return s
end

function from_assignment_match_iostream(lhs, rhs::ANY, linfo)
    return ""
end


text:
# test null hypothesis cases

using Test
using IncrementalInference
# going to introduce two new constraint types
# import IncrementalInference: getSample

@testset "Post 237 without nullhypo" begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)

addFactor!(fg, [:x0;], Prior(Normal()))
addFactor!(fg, [:x0;:x1], LinearRelative(Normal(10,1)))

solveTree!(fg)

pts_ = getBelief(fg, :x1) |> getPoints

@test 80 < sum(2 .< pts_ .< 18)

end

@testset "Post 237 nullhypo on prior" begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addFactor!(fg, [:x0;], Prior(Normal(10,1)), nullhypo=0.5)

solveTree!(fg)

pts_ = getBelief(fg, :x0) |> getPoints

@test 10 < sum(-8 .< pts_ .< 8) < 80
@test 30 < sum(5 .< pts_ .< 15) < 85

end

@testset "Post 237 with nullhypo test" begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)

addFactor!(fg, [:x0;], Prior(Normal()))
addFactor!(fg, [:x0;:x1], LinearRelative(Normal(10,1)), nullhypo=0.5)

pts = approxConv(fg, :x0x1f1, :x1)

@test 20 < sum(pts .< 5)
@test 20 < sum(5 .< pts .< 15)

solveTree!(fg)

pts2 = getBelief(fg, :x1) |> getPoints

@test 30 < sum(8 .< pts2 .< 12) <= 80
@test 80 < sum(-10 .< pts2 .< 30)

end

#


text:
# # [PDE Discovery : Heat equations](@id heat_equation)
#
# Similar to the discovery of ODE and DAE systems, sparse regression can be used to discover the underlying
# equations for partial differential equations as well as proposed by [PDEFind](https://www.science.org/doi/10.1126/sciadv.1602614). The following example shows how this
# can be achieved with DataDrivenDiffEq.jl, using the heat equation with Dirichlet boundary conditions
# with the analytical soultion $u(x,t) = sin(2\pi x) exp^{-(2\pi t)^2t}$.
# We start by defining the system and generate some data.

using DataDrivenDiffEq
using LinearAlgebra
using ModelingToolkit
using OrdinaryDiffEq
using DiffEqOperators
#md using Plots

u_analytic(x, t) = sin(2*π*x) * exp(-t*(2*π)^2)
nknots = 100
h = 1.0/(nknots+1)
knots = range(h, step=h, length=nknots)
ord_deriv = 2
ord_approx = 2

const bc = Dirichlet0BC(Float64)
const Δ = CenteredDifference(ord_deriv, ord_approx, h, nknots)

t0 = 0.0
t1 = 1.0
u0 = u_analytic.(knots, t0)

step(u,p,t) = Δ*bc*u
prob = ODEProblem(step, u0, (t0, t1))
alg = KenCarp4()
de_solution = solve(prob, alg);

#md plot(de_solution, legend = nothing)

# Using DiffEqOperators.jl, we can define the difference operators up to order $n =4$ and vectorize the result.

∂U = reduce(vcat, map(1:4) do n
    δ = CenteredDifference(n, ord_approx, h, nknots)
    reshape(
        δ*bc*Array(de_solution), 1, prod(size(de_solution))
        )
    end)

# Next we collect the discretized data samples, their time derivatives and define a DataDrivenProblem

U = reshape(Array(de_solution), 1, prod(size(de_solution)))
∂ₜU = reshape(Array(de_solution(de_solution.t, Val{1})), 1, prod(size(de_solution)))

problem = DataDrivenProblem(
    U, DX = ∂ₜU, U = ∂U
)

# We choose to model the spatial derivatives $\frac{d^n u}{dx^n}$ as exegenous signals (controls), which are directly substituted into the
# Basis defined in the following

@parameters t x
@variables u(x, t)
∂u = map(1:4) do n
    d = Differential(x)^n
    d(u)
end

basis = Basis([monomial_basis([u], 5);∂u], [u], independent_variable = t, controls = ∂u);
println(basis) #hide

# Afterwards, we define a sampler for the available data which performs a 80-20 train-test split and partions the training
# data into 10 batches and solve the sparse regression using STLSQ.

sampler = DataSampler(Split(ratio = 0.8), Batcher(n = 10))
solution = solve(problem, basis, STLSQ(1e-2:1e-2:5e-1), sampler = sampler, by = :best)

# As we can see, the heat equation is recovered correclty
result(solution)
println(result(solution)) #hide

#md # ## [Copy-Pasteable Code](@id cartpole_copy_paste)
#md #
#md # ```julia
#md # @__CODE__
#md # ```

## Test #src
for r_ in [solution] #src
    @test all(l2error(r_) .< 0.5) #src
    @test all(aic(r_) .> 1e3) #src
    @test all(determination(r_) .>= 0.9) #src
end #src


text:
################################################################################
### For generic methods
################################################################################

# Heuristic for selecting the level with the (presumably) most reliable
# standard error estimate:
# Take the highest lvl with at least 32 bins.
# (Chose 32 based on https://doi.org/10.1119/1.3247985)
function _reliable_level(B::LogBinner{T,N})::Int64 where {T,N}
    isempty(B) && (return 1)                # results in NaN in std_error
    i = findlast(x -> x.count >= 32, B.accumulators)
    return something(i, 1)
end

_eachlevel(B::LogBinner) = 1:findlast(acc -> acc.count > 1, B.accumulators)

################################################################################
### Statistics
################################################################################

"""
    mean(B::LogBinner[, lvl])

Calculates the mean for a given level in the Binning Analysis.
"""
mean(B::LogBinner, lvl = 1) = mean(B.accumulators[lvl])

"""
    var(B::LogBinner[, lvl])

Calculates the variance of a given level in the Binning Analysis.
"""
var(B::LogBinner, lvl = _reliable_level(B)) = var(B.accumulators[lvl])

"""
    varN(B::LogBinner[, lvl])

Calculates the variance/N of a given level in the Binning Analysis.
"""
varN(B::LogBinner, lvl = _reliable_level(B)) = varN(B.accumulators[lvl])

################################################################################
### Experimental
################################################################################

"""
    convergence(B::LogBinner, lvl)

Computes the difference between the variance of this lvl and the last,
normalized to the last lvl. If this value tends to 0, the Binning Analysis has
converged.
"""
function convergence(B::LogBinner) end

function convergence(B::LogBinner{T,N}, lvl = _reliable_level(B)) where {N, T <: Number}
    return abs((varN(B, lvl+1) - varN(B, lvl)) / varN(B, lvl))
end
function convergence(B::LogBinner{T,N}, lvl = _reliable_level(B)) where {N, T <: AbstractArray}
    return mean(abs.((varN(B, lvl+1) .- varN(B, lvl)) ./ varN(B, lvl)))
end

"""
    has_converged(B::LogBinner, lvl[, threshhold = 0.05])

Returns true if the Binning Analysis has converged for a given lvl.
"""
function has_converged(B::LogBinner, lvl = _reliable_level(B), threshhold::Float64 = 0.05)
    return convergence(B, lvl) <= threshhold
end

max_stars_repo_name: crstnbr/BinningAnalysis

max_issues_repo_name: crstnbr/BinningAnalysis

max_forks_repo_name: crstnbr/BinningAnalysis


text:
module InformativePrior

using Franklin

"""
Adapted from Franklin code at src/manager/post_processing.jl:126

This is a simple wrapper doing a git commit and git push without much
fanciness. It assumes the current directory is a git folder.
It also fixes all links if you specify `prepath` (or if it's set in
`config.md`).

**Keyword arguments**

* `prerender=true`: prerender javascript before pushing see
                     [`optimize`](@ref)
* `minify=true`:    minify output before pushing see [`optimize`](@ref)
* `nopass=false`:   set this to true if you have already run `optimize`
                     manually.
* `prepath=""`:     set this to something like "project-name" if it's a
                     project page
* `message="franklin-update"`: add commit message.
* `cleanup=true`:   whether to cleanup environment dictionaries (should
                     stay true).
* `final=nothing`:  a function `()->nothing` to execute last, before doing
                     the git push. It can be used to refresh a Lunr index,
                     generate notebook files with Literate, etc, ...
                     You might want to compose `fdf_*` functions that are
                     exported by Franklin (or imitate those). It can
                     access GLOBAL_VARS.
"""
function deploy(; prerender::Bool=true, minify::Bool=true, nopass::Bool=false,
                   prepath::String="", message::String="franklin-update",
                   cleanup::Bool=true, final::Union{Nothing,Function}=nothing,
                   do_push::Bool=true)::Nothing
    succ = true
    if !isempty(prepath) || !nopass
        # no cleanup so that can access global page variables in final step
        succ = optimize(prerender=prerender, minify=minify,
                        sig=true, prepath=prepath, cleanup=false)
    end
    if succ
        # final hook
        final === nothing || final()
        # --------------------------
        run(`rm -rf __production`)
        run(`cp -r __site/ __production`)

        # Push to git (publication)
        start = time()
        pubmsg = rpad("→ Pushing updates with git...", 35)
        print(pubmsg)
        try
            run(`git add -A `)
            wait(run(`git commit -m "$message" --quiet`; wait=false))
            if do_push
                run(`git push --quiet`)
            end
            print_final(pubmsg, start)
        catch e
            println("✘ Could not push updates, verify your connection " *
                    "and/or try manually.\n")
            @show e
        end
    else
        println("✘ Something went wrong in the optimisation step. " *
                "Not pushing updates.")
    end
    return nothing
end

end

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
const LONG_MAX = Base.typemin(Int32)
const AttrNumber = Int16
const FdwRoutine = Cvoid
const BlockNumber = UInt32
const HTAB = Cvoid
const OffsetNumber = UInt16
const MemoryContext = Ptr{Cvoid} # Ptr{MemoryContextData}
const PLpgSQL_function = Cvoid
const bool = Cuchar
const TupleTableSlot = Cvoid
const dsa_area = Cvoid
const dsa_pointer = Ptr{Cvoid} # UInt64
const ItemPointer = Ptr{Cvoid} # UInt64
const int8 = UInt8
const uint8 = UInt8
const int16 = UInt16
const uint16 = UInt16
const int32 = Int32
const uint32 = UInt32
const bits8 = uint8
const bits16 = uint16
const bits32 = uint32
const Size = UInt64
const Index = uint32
const Offset = int32
const UINT_MAX = Base.typemax(UInt64)
const Datum = UInt64
#=
union ValUnion
{
    long        ival;       /* machine integer */
    char       *str;        /* string */
}           val;
=#
const ValUnion = Ptr{Cvoid}

# struct ANONYMOUS1_data
#     #=
#     union
#     {
#         void       *ptr_value;
#         int         int_value;
#         Oid         oid_value;
#     }           data;
#     =#
#     ptr_value::Ptr{Cvoid} # Union{Ptr{List}, Cint, Oid=Cint}
# end


text:
"""
A = JopDiagonal(spc, d)
A = JopDiagonal(d)

where `spc::JetSpace` is the domain/range of `A`, and `d::Array` or `d::Number`
is the diagonal.  If `d<:Number`, then the diagonal of the matrix is constant
and one must specify `spc`.

# Examples
```julia
A = JopDiagonal(spc, d)
```
where `spc::JetSpace` is the domain and range of `A`, and `d<:Array` or
`d<:Number`

```julia
A = JotOpDiagonal(d)
```
where `d<:AbstractArray`.  The domain and range of `A` are determined by the
size and type of `d`.

```julia
A = JotOpDiagonal([1.0, 2.0, 3.0])
m = ones(domain(A))
d = A*m # d = [1.0 ; 2.0 ; 3.0]
```
"""
function JopDiagonal(diagonal::AbstractArray{T}) where {T}
    spc = space(diagonal)
    JopLn(df! = JopDiagonal_df!, df′! = JopDiagonal_df′!, dom = spc, rng = spc, s = (diagonal=diagonal,))
end
JopDiagonal(spc::JetAbstractSpace, diagonal) = JopLn(df! = JopDiagonal_df!, df′! = JopDiagonal_df′!, dom = spc, rng = spc, s = (diagonal=diagonal,))
export JopDiagonal

JopDiagonal_df!(d, m; diagonal, kwargs...) = d .= diagonal .* m
JopDiagonal_df′!(m, d; diagonal, kwargs...) = m .= conj(diagonal) .* d


text:
include("..\\src\\introducao.jl")

module IntroducaoTeste
  include("error_handler_pt_br.jl")
end

using Base.Test

import Introducao

Test.with_handler(custom_handler) do
  @test soma(1,2) == 2
end

max_stars_repo_name: raphaelpanta/julia-lang-exemplos

max_issues_repo_name: raphaelpanta/julia-lang-exemplos

max_forks_repo_name: raphaelpanta/julia-lang-exemplos


text:
using OrdinaryDiffEq, StochasticDiffEq, Flux, DiffEqFlux, Test

xs = Float32.(hcat([0.; 0.], [1.; 0.], [2.; 0.]))
tspan = (0.0f0,25.0f0)
dudt = Chain(Dense(2,50,tanh),Dense(50,2))

@testset "only end" begin
    only_end = neural_ode(dudt,xs,tspan,Tsit5(),save_everystep=false,save_start=false)
    @test_broken size(only_end)[end]==1
    CuArrays.allowscalar(false)
    only_end = neural_ode(dudt,xs,tspan,Tsit5(),save_everystep=false,save_start=false)
    @test_broken size(only_end)[end]==1
end


text:
module CellArrayReindex

using Gridap
using Gridap.CellValuesGallery
using Gridap.CachedArrays
using Base: @pure

import Gridap: reindex
import Base: iterate
import Base: length
import Base: size
import Base: getindex
import Base: IndexStyle

function reindex(values::IndexCellArray, indices::CellValue{<:IndexLike})
  IterCellArrayFromReindex(values,indices)
end

function reindex(values::IndexCellArray, indices::IndexCellValue{<:IndexLike})
  IndexCellArrayFromReindex(values,indices)
end

struct IterCellArrayFromReindex{T,N,V,I} <: IterCellValue{CachedArray{T,N,Array{T,N}}}
  values::V
  indices::I
end

function IterCellArrayFromReindex(
  values::IndexCellArray{T,N}, indices::CellValue{<:IndexLike}) where {T,N}
  V = typeof(values)
  I = typeof(indices)
  IterCellArrayFromReindex{T,N,V,I}(values,indices)
end

@inline function iterate(s::IterCellArrayFromReindex{T,N}) where {T,N}
  inext = iterate(s.indices)
  cache = CachedArray(T,N)
  _iterate(s,inext,cache)
end

@inline function iterate(s::IterCellArrayFromReindex,state)
  istate, cache = state
  inext = iterate(s.indices,istate)
  _iterate(s,inext,cache)
end

@inline function _iterate(s,inext,cache)
  if inext === nothing; return nothing end
  i, istate = inext
  v = s.values[i]
  _update_cache!(cache,v)
  state = (istate, cache)
  (cache, state)
end

function _update_cache!(cache,v)
  setsize!(cache,size(v))
  for (j,vi) in zip(eachindex(cache),v)
    cache[j] = vi
  end
end

length(s::IterCellArrayFromReindex) = length(s.indices)

struct IndexCellArrayFromReindex{T,N,D,V,I} <: IndexCellValue{CachedArray{T,N,Array{T,N}},D}
  values::V
  indices::I
  cache::CachedArray{T,N,Array{T,N}}
end

function IndexCellArrayFromReindex(
  values::IndexCellArray{T,N}, indices::IndexCellValue{<:IndexLike,D}) where {D,T,N}
  V = typeof(values)
  I = typeof(indices)
  cache = CachedArray(T,N)
  IndexCellArrayFromReindex{T,N,D,V,I}(values,indices,cache)
end

function getindex(
  s::IndexCellArrayFromReindex{T,D}, i::Vararg{<:Integer,D}) where {T,D}
  j = s.indices[i...]
  v = s.values[j]
  _update_cache!(s.cache,v)
  s.cache
end

size(s::IndexCellArrayFromReindex) = size(s.indices)

@pure function IndexStyle(
  ::Type{<:IndexCellArrayFromReindex{T,N,D,V,I}}) where {T,N,D,V,I}
  IndexStyle(I)
end

end # module


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using JuMP, LinearAlgebra, Geodesy, Cbc, Clp, ProgressBars, Printf, DataStructures

function get_solution(model::JuMP.Model; marginal_costs = true)
    graph, instance = model[:graph], model[:instance]
    T = instance.time

    output = OrderedDict(
        "Plants" => OrderedDict(),
        "Products" => OrderedDict(),
        "Costs" => OrderedDict(
            "Fixed operating (\$)" => zeros(T),
            "Variable operating (\$)" => zeros(T),
            "Opening (\$)" => zeros(T),
            "Transportation (\$)" => zeros(T),
            "Disposal (\$)" => zeros(T),
            "Expansion (\$)" => zeros(T),
            "Storage (\$)" => zeros(T),
            "Total (\$)" => zeros(T),
        ),
        "Energy" =>
            OrderedDict("Plants (GJ)" => zeros(T), "Transportation (GJ)" => zeros(T)),
        "Emissions" => OrderedDict(
            "Plants (tonne)" => OrderedDict(),
            "Transportation (tonne)" => OrderedDict(),
        ),
    )

    plant_to_process_node = OrderedDict(n.location => n for n in graph.process_nodes)
    plant_to_shipping_nodes = OrderedDict()
    for p in instance.plants
        plant_to_shipping_nodes[p] = []
        for a in plant_to_process_node[p].outgoing_arcs
            push!(plant_to_shipping_nodes[p], a.dest)
        end
    end

    # Products
    if marginal_costs
        for n in graph.collection_shipping_nodes
            location_dict = OrderedDict{Any,Any}(
                "Marginal cost (\$/tonne)" => [
                    round(abs(JuMP.shadow_price(model[:eq_balance][n, t])), digits = 2) for t = 1:T
                ],
                "Latitude (deg)" => n.location.latitude,
                "Longitude (deg)" => n.location.longitude,
                "Amount (tonne)" => n.location.amount,
            )
            if n.product.name ∉ keys(output["Products"])
                output["Products"][n.product.name] = OrderedDict()
            end
            output["Products"][n.product.name][n.location.name] = location_dict
        end
    end

    # Plants
    for plant in instance.plants
        skip_plant = true
        process_node = plant_to_process_node[plant]
        plant_dict = OrderedDict{Any,Any}(
            "Input" => OrderedDict(),
            "Output" =>
                OrderedDict("Send" => OrderedDict(), "Dispose" => OrderedDict()),
            "Input product" => plant.input.name,
            "Total input (tonne)" => [0.0 for t = 1:T],
            "Total output" => OrderedDict(),
            "Latitude (deg)" => plant.latitude,
            "Longitude (deg)" => plant.longitude,
            "Capacity (tonne)" =>
                [JuMP.value(model[:capacity][process_node, t]) for t = 1:T],
            "Opening cost (\$)" => [
                JuMP.value(model[:open_plant][process_node, t]) *
                plant.sizes[1].opening_cost[t] for t = 1:T
            ],
            "Fixed operating cost (\$)" => [
                JuMP.value(model[:is_open][process_node, t]) *
                plant.sizes[1].fixed_operating_cost[t] +
                JuMP.value(model[:expansion][process_node, t]) *
                slope_fix_oper_cost(plant, t) for t = 1:T
            ],
            "Expansion cost (\$)" => [
                (
                    if t == 1
                        slope_open(plant, t) * JuMP.value(model[:expansion][process_node, t])
                    else
                        slope_open(plant, t) * (
                            JuMP.value(model[:expansion][process_node, t]) -
                            JuMP.value(model[:expansion][process_node, t-1])
                        )
                    end
                ) for t = 1:T
            ],
            "Process (tonne)" =>
                [JuMP.value(model[:process][process_node, t]) for t = 1:T],
            "Variable operating cost (\$)" => [
                JuMP.value(model[:process][process_node, t]) *
                plant.sizes[1].variable_operating_cost[t] for t = 1:T
            ],
            "Storage (tonne)" =>
                [JuMP.value(model[:store][process_node, t]) for t = 1:T],
            "Storage cost (\$)" => [
                JuMP.value(model[:store][process_node, t]) * plant.storage_cost[t]
                for t = 1:T
            ],
        )
        output["Costs"]["Fixed operating (\$)"] += plant_dict["Fixed operating cost (\$)"]
        output["Costs"]["Variable operating (\$)"] +=
            plant_dict["Variable operating cost (\$)"]
        output["Costs"]["Opening (\$)"] += plant_dict["Opening cost (\$)"]
        output["Costs"]["Expansion (\$)"] += plant_dict["Expansion cost (\$)"]
        output["Costs"]["Storage (\$)"] += plant_dict["Storage cost (\$)"]

        # Inputs
        for a in process_node.incoming_arcs
            vals = [JuMP.value(model[:flow][a, t]) for t = 1:T]
            if sum(vals) <= 1e-3
                continue
            end
            skip_plant = false
            dict = OrderedDict{Any,Any}(
                "Amount (tonne)" => vals,
                "Distance (km)" => a.values["distance"],
                "Latitude (deg)" => a.source.location.latitude,
                "Longitude (deg)" => a.source.location.longitude,
                "Transportation cost (\$)" =>
                    a.source.product.transportation_cost .* vals .* a.values["distance"],
                "Transportation energy (J)" =>
                    vals .* a.values["distance"] .* a.source.product.transportation_energy,
                "Emissions (tonne)" => OrderedDict(),
            )
            emissions_dict = output["Emissions"]["Transportation (tonne)"]
            for (em_name, em_values) in a.source.product.transportation_emissions
                dict["Emissions (tonne)"][em_name] =
                    em_values .* dict["Amount (tonne)"] .* a.values["distance"]
                if em_name ∉ keys(emissions_dict)
                    emissions_dict[em_name] = zeros(T)
                end
                emissions_dict[em_name] += dict["Emissions (tonne)"][em_name]
            end
            if a.source.location isa CollectionCenter
                plant_name = "Origin"
                location_name = a.source.location.name
            else
                plant_name = a.source.location.plant_name
                location_name = a.source.location.location_name
            end

            if plant_name ∉ keys(plant_dict["Input"])
                plant_dict["Input"][plant_name] = OrderedDict()
            end
            plant_dict["Input"][plant_name][location_name] = dict
            plant_dict["Total input (tonne)"] += vals
            output["Costs"]["Transportation (\$)"] += dict["Transportation cost (\$)"]
            output["Energy"]["Transportation (GJ)"] +=
                dict["Transportation energy (J)"] / 1e9
        end

        plant_dict["Energy (GJ)"] = plant_dict["Total input (tonne)"] .* plant.energy
        output["Energy"]["Plants (GJ)"] += plant_dict["Energy (GJ)"]

        plant_dict["Emissions (tonne)"] = OrderedDict()
        emissions_dict = output["Emissions"]["Plants (tonne)"]
        for (em_name, em_values) in plant.emissions
            plant_dict["Emissions (tonne)"][em_name] =
                em_values .* plant_dict["Total input (tonne)"]
            if em_name ∉ keys(emissions_dict)
                emissions_dict[em_name] = zeros(T)
            end
            emissions_dict[em_name] += plant_dict["Emissions (tonne)"][em_name]
        end

        # Outputs
        for shipping_node in plant_to_shipping_nodes[plant]
            product_name = shipping_node.product.name
            plant_dict["Total output"][product_name] = zeros(T)
            plant_dict["Output"]["Send"][product_name] = product_dict = OrderedDict()

            disposal_amount = [JuMP.value(model[:dispose][shipping_node, t]) for t = 1:T]
            if sum(disposal_amount) > 1e-5
                skip_plant = false
                plant_dict["Output"]["Dispose"][product_name] =
                    disposal_dict = OrderedDict()
                disposal_dict["Amount (tonne)"] =
                    [JuMP.value(model[:dispose][shipping_node, t]) for t = 1:T]
                disposal_dict["Cost (\$)"] = [
                    disposal_dict["Amount (tonne)"][t] *
                    plant.disposal_cost[shipping_node.product][t] for t = 1:T
                ]
                plant_dict["Total output"][product_name] += disposal_amount
                output["Costs"]["Disposal (\$)"] += disposal_dict["Cost (\$)"]
            end

            for a in shipping_node.outgoing_arcs
                vals = [JuMP.value(model[:flow][a, t]) for t = 1:T]
                if sum(vals) <= 1e-3
                    continue
                end
                skip_plant = false
                dict = OrderedDict(
                    "Amount (tonne)" => vals,
                    "Distance (km)" => a.values["distance"],
                    "Latitude (deg)" => a.dest.location.latitude,
                    "Longitude (deg)" => a.dest.location.longitude,
                )
                if a.dest.location.plant_name ∉ keys(product_dict)
                    product_dict[a.dest.location.plant_name] = OrderedDict()
                end
                product_dict[a.dest.location.plant_name][a.dest.location.location_name] =
                    dict
                plant_dict["Total output"][product_name] += vals
            end
        end

        if !skip_plant
            if plant.plant_name ∉ keys(output["Plants"])
                output["Plants"][plant.plant_name] = OrderedDict()
            end
            output["Plants"][plant.plant_name][plant.location_name] = plant_dict
        end
    end

    output["Costs"]["Total (\$)"] = sum(values(output["Costs"]))
    return output
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
module ImputedGenotypes
	using BGZFStreams

	export r², likelihoods_to_dosage, mean_gp, best_gp, confidence_weighted_mean_gp, merge_imputed

	typealias GenotypeProbabilities Array{Tuple{Float64,Float64,Float64},1}
	typealias GenotypeProbability Tuple{Float64,Float64,Float64}

	probs(pr::Array{Tuple{Float64,Float64,Float64},1}) = normalize([foldr((x,y)->(x[1]+y[1],x[2]+y[2], x[3]+y[3]), (.0,.0,.0), pr)...],1)
	freqA(g) = 0.5*probs(g)[2] + probs(g)[1]

	function r²(ps::Array{Tuple{Float64,Float64,Float64},1})
		n = length(ps)

		esum = e2sum = fsum = 0.0

		for l in ps
			norm = sum(l)
			l = [ i/norm for i in l ]
			esum += l[2] + 2l[3]
			e2sum += (l[2] + 2l[3])^2
			fsum += l[2] + 4l[3]
		end

		θ = esum / 2n

		(1.0 > freqA(ps) > 0.0) ?  abs(1 - (fsum - e2sum) / (2n * θ * (1.0 - θ))) : 1.0
	end

	likelihoods_to_dosage(p_aa, p_ab, p_bb) = 2p_bb + p_ab
	likelihoods_to_dosage{T<:AbstractFloat}(triple::Tuple{T,T,T}) = likelihoods_to_dosage(triple[1],triple[2],triple[3])

	function max_index(p)
		max_index = 0
		max_probability = 0.0
		for i in 1:length(p)
			if p[i] > max_probability
				max_index = i
				max_probability = p[i]
			end
		end
		max_index
	end

	max_prob(p) = p[max_index(p)]

	putf(f) = @sprintf("%0.3f",f)

	alt_allele_count(p::Tuple{Float64,Float64,Float64}) = max_index(p)-1
	alt_allele_count(p::Array{Tuple{Float64,Float64,Float64},1}) = sum(map(alt_allele_count,p))

	GT(p::Tuple{Float64,Float64,Float64}) = ("0/0","0/1","1/1")[max_index(p)]
	DS(p::Tuple{Float64,Float64,Float64}) = putf(likelihoods_to_dosage(p[1],p[2],p[3]))
	GP(p::Tuple{Float64,Float64,Float64}) = join(map(putf,p),',')

	INFO(l)=string("AC=",length(l)*2,";AN=",alt_allele_count(l),";R2=", r²(l))

	best_gp(ps) = ps[max_index(map(max_prob,ps))]

	mean_gp(ps) = tuple(probs(ps)...)

	function confidence_weighted_mean_gp(ps)
		cs = map(p -> foldl(max,p) - foldl(min,p), ps)
		tuple(probs([ (cs[i]*ps[i][1], cs[i]*ps[i][2], cs[i]*ps[i][3]) for i in 1:length(ps) ])...)
	end

	best_info(ps) = () # Only a Placeholder

	function merge_snp(fields, mergefun)
		@assert length(fields) > 0

		println(join(fields[1][1:4],' '))

		gpidx = [ first(find(x -> x=="GP"||x=="GL", split(f[9],':'))) for f in fields ]

		# Create an array Array{Tuple{Float64,Float64,Float64},1} for each individual
		best_likelihoods = []
		let snp_probs(i,indv) = (map(x->parse(Float64,x),split(split(fields[i][indv],':')[gpidx[i]],','))...)
			# Not very multi-dispatchy, but ...
			if mergefun == best_info
				# This needs access to others genotypes as well
				# likelihoods is an array  of Array{Tuple{Float64,Float64,Float64},1} for each method
				likelihoods = [ [ snp_probs(idx,indv) for indv in 10:length(first(fields)) ] for idx in 1:length(fields) ]
				# Select likehoods of method with highest overall INFO score
				best_likelihoods = likelihoods[max_index(map(r², likelihoods))]
			else
				# Create an array Array{Tuple{Float64,Float64,Float64},1} for each individual:
				best_likelihoods = @parallel vcat for indv in 10:length(fields[1])
					# Create an array Array{Tuple{Float64,Float64,Float64},1} for each method and collapse with "mergefun"
					mergefun( filter(isna,[ snp_probs(i,indv) for i in 1:length(fields) ]) )
				end
			end
		end

		join(foldl(vcat,
			[fields[1][1:7],
			INFO(best_likelihoods),
			["GT:DS:GP"],
			map(p->join([GT(p),DS(p),GP(p)],':'), best_likelihoods),
			['\n']]),'\t')
	end

	vcf_header="""##fileformat=VCFv4.1
	##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
	##FORMAT=<ID=DS,Number=1,Type=Float,Description="Estimated Alternate Allele Dosage : [P(0/1)+2*P(1/1)]">
	##FORMAT=<ID=GP,Number=3,Type=Float,Description="Estimated Posterior Probabilities for Genotypes 0/0, 0/1 and 1/1 ">
	##INFO=<ID=R2,Number=1,Type=Float,Description="Estimated Imputation Accuracy">
	"""

	# Two vcf lines are considered to represent the same variant if
	# chromosome AND position AND ref+alt allele are the same (but order of alt/ref may be opposite)
	function vcfline_same_variant(f1,f2)
		f1[1] == f2[1] && f1[2] == f2[2] && sort(f1[3:4]) == sort(f2[3:4])
	end

	# Greater-than operator for vcf lines (split into fields)
	function vcfline_before(f1,f2)
		if f1[1] < f2[1] # Chromosome smaller
			true
		elseif f1[1] == f2[1] && f1[2] < f2[2] # Position smaller
			true
		elseif f1[1] == f2[1] && f1[2] == f2[2] && string(sort(f1[3:4])) < string(sort(f2[3:4])) # allele ordering
			true
		else
			false
		end
	end

	revidx(lines) = [ Dict(f[i] => i for i in 10:length(f)) for f in in [ split(line) for line in lines ] ]

	function merge_imputed(vcfs_files, out_file, mergefun)
		vcfs = [ BGZFStream(vcf) for vcf in vcf_files ]
		out = BGZFStream(out_file)

		lines = []
		states = []
		current_lines = []
		for vcf in vcfs
			lines_vcf = eachline(vcf)
			push!(lines,lines_vcf)
			push!(states,start(lines_vcf))
			push!(current_lines,"##")
		end

		# Read past header lines
		for i in 1:length(lines)
			while ismatch(r"^##", current_lines[i])
				current_lines[i], states[i] = next(lines[i], states[i])
			end
		end

		@assert (all(map(line->ismatch(r"^#CHROM", line),current_lines)))

		# Build sample indices
		sample_idx=revidx(current_lines)

		# TODO:
		# And idea for a strategy that requires minimal changes to merge_snp function
		# would be to "reconstruct" full lines with all samples in a paricular order
		# where some samples may be NA.
		# 2) As an optimization we may only have to do this when not all files contain
		# all samples in the same order...

		samples_in_order = sort(foldl(union,map(keys,sample_idx)))

		write(out,vcf_header)
		write(out,current_lines[1])

		for i in 1:length(lines)
			current_lines[i], states[i] = next(lines[i], states[i])
		end
		current_fields = [ split(line) for line in current_lines ]

		while all(map(x->!done(x,nothing), lines))
			# Find smallest position
			smallest =  1
			for i in 1:length(lines)
				if vcfline_before(current_fields[i], current_fields[smallest])
					smallest=i
				end
			end

			# Select only lines with smallest position for merging in this step and read next line from those files
			process_lines = []
			for i in 1:length(lines)
				if vcfline_same_variant(current_fields[i], current_fields[smallest])
					current_lines[i], states[i] = next(lines[i], states[i])
					current_fields[i] = split(current_lines[i])
					push!((i,process_lines), current_fields[i])
				end
			end

			# I am going to need the index for each process_lines (To reverse map individuals
			# I need to know hwich file they are from

			# TODO:
			# And idea for a strategy that requires minimal changes to merge_snp function
			# would be to "reconstruct" full lines with all samples in a paricular order
			# where some samples may be NA.

			proclines = []
			for (idx,line) in process_lines
				aline=[]
				sidx = sample_idx[idx]
				for i in 1:length(procline)
					if < 10 # Copy over non-genotype fields
						push!(aline, line[i])
					elseif haskey(samples[idx],samples_in_order[i-10])
						push!(aline,line[sample_idx[idx][samples_in_order[i-10]]])
					else
						push!(aline,NA)
					end
				end
				push!(proclines,aline)
			end

			write(out,merge_snp(proclines,mergefun))
		end

		# Close open files
		for vcf in vcfs
			close(vcf)
		end
		close(out)
	end

end

max_stars_repo_name: cth/merge-imputations

max_issues_repo_name: cth/merge-imputations

max_forks_repo_name: cth/merge-imputations


text:
export Response,
       start,
       flush,
       done,
       send,
       sentcl,
       headerssent,
       clean,
       getclean,
       data,
       buffered,
       chunk,
       chuncked,
       file

defaultresponseheaders() = headers([
    ("Server"            , "Julia/$(Base.VERSION) at $NAME/$VERSION"),
    ("Content-Type"      , "text/html; charset=utf-8"),
    ("Date"              , RFC1123_datetime())
])

type ResponseException <: Exception
    info
end

ResponseException(key::Union(STR_TYPE, (STR_TYPE...)), message::STR_TYPE; kwargs...) = begin
    all = append!(Any[(:key, key), (:message, message)], kwargs)
    ResponseException(Dict{Any, Any}(all))
end

type Response
    io::AbstractIOSocket
    status::Integer
    headers::Headers
    data
    finished::Bool
    protocol::Protocol
    headers_sent::Bool
    writer::Union(N_TYPE,Function)
    buffer::Union(N_TYPE, IO)
    content_length::Integer
    sent_bytes::Integer
    headers_size::Integer
end

const EMPTY_CL = -1

function startwriter(r::Response, data = N)
    if data != N && !chunked(r) && r.content_length <= 0 && r.buffer == N
        r.content_length = sizeof(data)
    end

    start(r)

    if data != N
        Base.write(r, data)
    end
end

Response(io::AbstractIOSocket, protocol::Protocol) = Response(
    io,
    protocol.default_status.default,
    defaultresponseheaders(),
    N,
    false,
    protocol,
    false,
    startwriter,
    N,
    EMPTY_CL,
    0,
    0
)

chunked(r::Response) = isa(r.io, Chunk.IOChunked)

function chunk(r::Response)
    if chunked(r)
        throw(ResponseException("chunk_enabled",
          "The Response.chunk has be previously enabled."))
    end
    r.io = Chunk.IOChunked(r.io.io; original=r.io)
    r
end

function sendheaders(r::Response)
    if r.headers_sent
        throw(ResponseException("headers_sent", "headers has be sent"))
    end

    r.headers_sent = true

    h_chunk, h_chunk_value = r.protocol.header_msg_chunk

    if r.content_length != EMPTY_CL
        r.headers[r.protocol.header_msg_size] = str(r.content_length)
    end

    if get(r.headers, h_chunk, "") == h_chunk_value && !chunked(r)
        chunk(r)
    end

    if chunked(r)
        r.headers[h_chunk] = h_chunk_value
    end

    if !(r.status in r.protocol.status_wo_msg) && !chunked(r) && r.content_length == EMPTY_CL
        r.headers[r.protocol.header_msg_size] = "0"
    end

    if haskey(r.headers, r.protocol.header_msg_size)
        r.content_length = integer(r.headers[r.protocol.header_msg_size])
    end

    r.protocol.prepare_res_headers(r)

    writedata(r, string(r.protocol.name, "/", r.protocol.version,
        " ", r.status, " ",
        r.protocol.status_codes[r.status], CRLF))

    for (header, value) in r.headers
        if isa(value, STR_TYPE)
            writedata(r, string(header, ": ", value, CRLF))
        else
            for v in value
                writedata(r, string(header, ": ", v, CRLF))
            end
        end
    end

    writedata(r, CRLF)
    r.headers_size = r.sent_bytes
end

bufferedwriter(r, data) = write(r.buffer, data)

function start(r::Response)
    if r.buffer == N
        sendheaders(r)
        r.writer = chunked(r) ? Chunk.write : (
            r.content_length > 0 ? writedatacheck: writedata
        )
    else
        r.writer = bufferedwriter
    end
end

Base.write(r::Response, data::STR_TYPE) = r.writer(r, data)
Base.write(r::Response, data::Array{Uint8}) = r.writer(r, data)

<<(r::Response, data::STR_TYPE) = r.writer(r, data)
<<(r::Response, data::Array{Uint8}) = r.writer(r, data)

>>(data::STR_TYPE, r::Response) = r.writer(r, data)
>>(data::Array{Uint8}, r::Response) = r.writer(r, data)

Base.flush(r::Response) = begin
    b = takebuf_array(r.buffer)

    if !r.headers_sent
        if !chunked(r) && r.content_length == EMPTY_CL
            r.content_length = sizeof(b)
        end

        sendheaders(r)

        if chunked(r)
            Chunk.write(r, b)
        else
            writedata(r, b)
        end
    elseif chunked(r)
        Chunk.write(r, b)
    else
        writedatacheck(r, b)
    end
end

function done(r::Response)
    if !r.headers_sent
        startwriter(r)
    end

    if r.buffer != N
        flush(r)
    end

    if r.data != N
        data = r.data
        r.data = N

        if isa(data, Union(STR_TYPE, Array{Uint8}))
            write(r, data)
        elseif isiter(data)
            for d in data
                write(r, d)
            end
        else
            throw(ResponseException("data_field", "Invalid Response.data
                type '$(typeof(data))'"), data=data)
        end
    end

    if chunked(r)
        r.sent_bytes += Chunk.done(r.io)
    end
end

clean(r::Response) = (buffer = IOBuffer())
getclean(r::Response) = takebufarray(r.buffer)
data(r::Response) = r.buffer.data
writedata(r::Response, data) = (s = write(r.io.io, data); r.sent_bytes += s; s)
sentcl(r::Response) = r.sent_bytes - r.headers_size

function writedatacheck(r::Response, data)
    left = r.content_length - sentcl(r)
    size = sizeof(data)

    if left - size < 0
        if left <= StreamReader.DEFAULT_PART_SIZE
            writedata(r, repeat("\0", left))
        else
            # the amount remaining to complete the content size is too large.
            # Let's break it small data blocks.
            # the num_parts > 1, because left > StreamReader.DEFAULT_PART_SIZE
            num_parts, parts_size, last_part_size = StreamReader.calculate(left)
            part_data = repeat("\0", parts_size)

            for i in 1:(num_parts-1) # skip last part
                writedata(r, part_data)
            end

            if last_part_size != parts_size
                writedata(r, part_data[1:last_part_size])
            else
                writedata(r, part_data)
            end
        end

        exceeded = size - left

        # args is (cl, size, exceeded)
        throw(ResponseException("exceeded_content_data",
            "the data size exceeded by $exceeded the total size of the content.";
            content_length=r.content_length, data_size=size, exceeded=exceeded,
            data=data))
    end
    writedata(r, data)
end

function sender(f::Function, r::Response, chunk::Bool = false)
    if chunk
        chunked(r) = true
    end

    start(r)
    f()
    done(r)
end

function buffered(r::Response)
    if r.buffer == N
        r.buffer = IOBuffer()
    end
end

show(io::IO, r::Response) = print(
    io,
    "Response(",
    r.status,
    " ",
    STATUS_CODES[r.status],
    ", ",
    length(r.headers),
    " Headers, ",
    isa(r.data, STR_TYPE) ? "$(sizeof(r.data)), Bytes in Body" : typeof(r.data),
    ")"
)

header(r::Response, key::STR_TYPE, value::STR_TYPE) = header(r.headers, key, value)

"""
Download file on Respose

flags:
  :fdw - Force File Download
"""
function file(r::Response, filename::STR_TYPE, bsize::Integer=0; download::Bool=false, dfilename::STR_TYPE="", kwargs...)
    finfo = fileinfo(filename)

    if finfo != N
        io = open(filename)
        ext, mime, fsize = finfo
        r.headers["Content-Type"] = mime

        if download
            if isempty(dfilename)
                dfilename = basename(filename)
            end
            r.headers["Content-Disposition"] = "attachment; filename=\"$dfilename\"";
        end

        if !chunked(r)
            r.content_length = fsize
        end

        r.data = StreamReader.IOReaderIterator(io, fsize, bsize; kwargs...)
    else
        r.status = r.protocol.default_status.not_found
        r.data = "Not Found - file $filename could not be found"
    end
end

Chunk.write(r::Response, data) = (r.sent_bytes += Chunk.write(r.io, data))


text:
# Some utility functions for extracting results.

type Sample
  weight :: Float64 # particle weight
  value :: Dict{Symbol,Any}
end

Base.getindex(s::Sample, v::Symbol) = Base.getindex(s.value, v)

type Chain
  weight :: Float64 # log model evidence
  value :: Array{Sample}
end
Chain() = Chain(0, Vector{Sample}())

function Base.getindex(c::Chain, v::Symbol)
  # This strange implementation is mostly to keep backward compatability.
  #  Needs some refactoring a better format for storing results is available.
  if v == :logevidence
    log(c.weight)
  else
    map((s)->Base.getindex(s, v), c.value)
  end
end

Base.push!(c::Chain, s::Sample) = push!(c.value, s)
# compute mean(f(x).w), where (x, w) is a weighted sample
Base.mean(c::Chain, v::Symbol, f::Function) = mapreduce((s)->f(s[v]).*s.weight, +, c.value)

Base.keys(p :: Particle) = keys(p.task.storage[:turing_predicts])
Base.values(p :: Particle) = values(p.task.storage[:turing_predicts])
Base.getindex(p :: Particle, args...) = getindex(p.task.storage[:turing_predicts], args...)

# ParticleContainer: particles ==> (weight, results)
function getsample(pc :: ParticleContainer, i :: Int, w :: Float64 = 0.)
  p = pc.vals[i]

  predicts = Dict{Symbol, Any}()
  for k in keys(p)
    predicts[k] = p[k]
  end
  return Sample(w, predicts)
end

function Chain(pc :: ParticleContainer)
  w = pc.logE
  chain = Array{Sample}(length(pc))
  Ws, z = weights(pc)
  s = map((i)->chain[i] = getsample(pc, i, Ws[i]), 1:length(pc))

  Chain(exp(w), s)
end

# tests
# tr = Turing.sampler.particles[1]
# tr = Chain(Turing.sampler.particles)


text:
# Render the widgets from WidgetsBase!

function jsrender(button::Button)
    return DOM.input(
        type = "button",
        value = button.content,
        onclick = js"update_obs($(button.value), true);";
        button.attributes...
    )
end

function jsrender(tf::TextField)
    return DOM.input(
        type = "textfield",
        value = tf.value,
        onchange = js"update_obs($(tf.value),  this.value);";
        tf.attributes...
    )
end

function jsrender(ni::NumberInput)
    return DOM.input(
        type = "number",
        value = ni.value,
        onchange = js"update_obs($(ni.value), parseFloat(this.value));";
        ni.attributes...
    )
end

function jsrender(slider::Slider)
    return DOM.input(
        type = "range",
        min = map(first, slider.range),
        max = map(last, slider.range),
        value = slider.value,
        step = map(step, slider.range),
        oninput = js"update_obs($(slider.value), parseFloat(value))";
        slider.attributes...
    )
end

const noUiSlider = Dependency(
    :noUiSlider,
    [
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.js",
        "https://cdn.jsdelivr.net/gh/leongersen/noUiSlider/distribute/nouislider.min.css"
    ]
)

function jsrender(session::Session, slider::RangeSlider)
    args = (slider.range, slider.connect, slider.orientation,
            slider.tooltips, slider.ticks)
    style = map(args...) do range, connect, orientation, tooltips, ticks
        value = slider.value[]
        return Dict(
            "range" => Dict(
                "min" => fill(range[1], length(value)),
                "max" => fill(range[end], length(value)),
            ),
            "step" => step(range),
            "start" => value,
            "connect" => connect,
            "tooltips" => tooltips,
            "direction" => "ltr",
            "orientation" => string(orientation),
            "pips" => isempty(ticks) ? nothing : ticks
        )
    end
    rangediv = DOM.div()
    create_slider = js"""function create_slider(style){
        var range = $(rangediv);
        range.noUiSlider.updateOptions(deserialize_js(style), true);
    }"""
    onload(session, rangediv, js"""function onload(range){
        var style = $(style[]);
        $(noUiSlider).create(range, style);
        range.noUiSlider.on('update', function (values, handle, unencoded, tap, positions){
            update_obs($(slider.value), [parseFloat(values[0]), parseFloat(values[1])]);
        });
    }""")
    onjs(session, style, create_slider)
    return rangediv
end

function JSServe.jsrender(tb::Checkbox)
    return DOM.input(
        type = "checkbox",
        checked = tb.value,
        onchange = js"update_obs($(tb.value), this.checked);";
        tb.attributes...
    )
end


text:
using SolidStateDetectors

sim = Simulation(SSD_examples[:InvertedCoax])

@info sim.detector.semiconductors[1].charge_density_model
@info supertype(typeof(sim.detector.semiconductors[1].charge_density_model))
# Custom made impurity profile:

struct CustomImpurityDensity{T} <: SolidStateDetectors.SolidStateDetectors.AbstractChargeDensityModel{T}
    # here you can add parameters of your profile (optional).
    # E.g.:
    offset::T # in 1/m^3
end

custom_impurity_density = CustomImpurityDensity{Float32}(1e16)

# overwrite the model defined in the config file
sim.detector.semiconductors[1].charge_density_model = custom_impurity_density

# calculate_electric_potential!(sim) will throw an error that a function is not defined
# since the packages needs to know how to evaluate the density at a certain point for the `CustomImpurityDensity`

function SolidStateDetectors.get_charge_density(cid::CustomImpurityDensity{T}, pt::CylindricalPoint{T}) where {T}
    # The return value must be in SI units: 1/m^3
    # dummy profile
    z = pt[3] # pt is in SI units. pt[1]: m, pt[2] = rad, pt[3] = m
    if z > 0.001 # > 1mm
        return cid.offset
    else
        0
    end
end

calculate_electric_potential!(sim)

# In case you are working in cartesian coordinates you will need to define it for cartesian points:
function SolidStateDetectors.get_charge_density(cid::CustomImpurityDensity{T}, pt::CartesianPoint{T}) where {T}
    # The return value must be in SI units: 1/m^3
    # dummy profile
    z = pt[3] # pt is in SI units. pt[1]: m, pt[2] = rad, pt[3] = m
    if z > 0.001 # > 1mm
        return cid.offset
    else
        0
    end
end

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
import Base: <=

typealias Interval{T} Tuple{T,T}

"""
```
Transform
```

Subtypes of the abstract Transform type indicate how a `Parameter`'s
value is transformed from model space (which can be bounded to a
limited section of the real line) to the entire real line (which is
necessary for mode-finding using csminwel). The transformation is
performed by the `transform_to_real_line` function, and is reversed by the
`transform_to_model_space` function.
"""
abstract Transform

immutable Untransformed <: Transform end
immutable SquareRoot    <: Transform end
immutable Exponential   <: Transform end

Base.show(io::IO, t::Untransformed) = @printf io "x -> x\n"
Base.show(io::IO, t::SquareRoot)    = @printf io "x -> (a+b)/2 + (b-a)/2*c*x/sqrt(1 + c^2 * x^2)\n"
Base.show(io::IO, t::Exponential)   = @printf io "x -> b + (1/c) * log(x-a)\n"

Base.call(t::Untransformed) = Untransformed()
Base.call(t::SquareRoot) = SquareRoot()
Base.call(t::Exponential) = Exponential()

"""
```
AbstractParameter{T<:Number}
```

The AbstractParameter type is the common supertype of all model
parameters, including steady-state values.  Its subtype structure is
as follows:

-`AbstractParameter{T<:Number}`: The common abstract supertype for all parameters.
    -`Parameter{T<:Number, U<:Transform}`: The abstract supertype for parameters that are directly estimated.
        -`UnscaledParameter{T<:Number, U:<Transform}`: Concrete type for parameters that do not need to be scaled for equilibrium conditions.
        -`ScaledParameter{T<:Number, U:<Transform}`: Concrete type for parameters that are scaled for equilibrium conditions.
    -`SteadyStateParameter{T<:Number}`: Concrete type for steady-state parameters.
"""
abstract AbstractParameter{T<:Number}

"""
```
Parameter{T<:Number, U<:Transform} <: AbstractParameter{T}
```

The Parameter type is the common supertype of time-invariant, non-steady-state model
parameters. It has 2 subtypes, `UnscaledParameter` and `ScaledParameter`.
`ScaledParameter`s are parameters whose values are scaled when used in the model's
equilibrium conditions. The scaled value is stored for convenience, and udpated when the
parameter's value is updated.
"""
abstract Parameter{T,U<:Transform} <: AbstractParameter{T}

typealias ParameterVector{T} Vector{AbstractParameter{T}}
typealias NullablePrior      Nullable{ContinuousUnivariateDistribution}

"""
```
UnscaledParameter{T<:Number,U<:Transform} <: Parameter{T,U}
```

Time-invariant model parameter whose value is used as-is in the model's equilibrium
conditions.

#### Fields
- `key::Symbol`: Parameter name. For maximum clarity, `key`
  should conform to the guidelines established in the DSGE Style Guide.
- `value::T`: Parameter value. Initialized in model space (guaranteed
  to be between `valuebounds`), but can be transformed between model
  space and the real line via calls to `transform_to_real_line` and
`transform_to_model_space`.
- `valuebounds::Interval{T}`: Bounds for the parameter's value in model space.
- `transform_parameterization::Interval{T}`: Parameters used to
  transform `value` between model space and the real line.
- `transform::U`: Transformation used to transform `value` between
  model space and real line.
- `prior::NullablePrior`: Prior distribution for parameter value.
- `fixed::Bool`: Indicates whether the parameter's value is fixed rather than estimated.
- `description::AbstractString`:  A short description of the parameter's economic
  significance.
- `tex_label::AbstractString`: String for printing the parameter name to LaTeX.
"""
type UnscaledParameter{T,U} <: Parameter{T,U}
    key::Symbol
    value::T                                # parameter value in model space
    valuebounds::Interval{T}                # bounds of parameter value
    transform_parameterization::Interval{T} # parameters for transformation
    transform::U                            # transformation between model space and real line for optimization
    prior::NullablePrior                    # prior distribution
    fixed::Bool                             # is this parameter fixed at some value?
    description::AbstractString
    tex_label::AbstractString               # LaTeX label for printing
end

"""
```
ScaledParameter{T,U} <: Parameter{T,U}
```

Time-invariant model parameter whose value is scaled for use in the model's equilibrium
conditions.

#### Fields

- `key::Symbol`: Parameter name. For maximum clarity, `key`
  should conform to the guidelines established in the DSGE Style Guide.
- `value::T`: The parameter's unscaled value. Initialized in model
  space (guaranteed to be between `valuebounds`), but can be
  transformed between model space and the real line via calls to
  `transform_to_real_line` and `transform_to_model_space`.
- `scaledvalue::T`: Parameter value scaled for use in `eqcond.jl`
- `valuebounds::Interval{T}`: Bounds for the parameter's value in model space.
- `transform_parameterization::Interval{T}`: Parameters used to
  transform `value` between model space and the real line.
- `transform::U`: The transformation used to convert `value` between model space and the
  real line, for use in optimization.
- `prior::NullablePrior`: Prior distribution for parameter value.
- `fixed::Bool`: Indicates whether the parameter's value is fixed rather than estimated.
- `scaling::Function`: Function used to scale parameter value for use in equilibrium
  conditions.
- `description::AbstractString`: A short description of the parameter's economic
  significance.
- `tex_label::AbstractString`: String for printing parameter name to LaTeX.
"""
type ScaledParameter{T,U} <: Parameter{T,U}
    key::Symbol
    value::T
    scaledvalue::T
    valuebounds::Interval{T}
    transform_parameterization::Interval{T}
    transform::U
    prior::NullablePrior
    fixed::Bool
    scaling::Function
    description::AbstractString
    tex_label::AbstractString
end

"""
```
SteadyStateParameter{T} <: AbstractParameter{T}
```

Steady-state model parameter whose value depends upon the value of other (non-steady-state)
`Parameter`s. `SteadyStateParameter`s must be constructed and added to an instance of a
model object `m` _after_ all other model `Parameter`s have been defined. Once added to `m`,
`SteadyStateParameter`s are stored in `m.steady_state`. Their values are calculated and set
by `steadystate!(m)`, rather than being estimated directly. `SteadyStateParameter`s do not
require transformations from the model space to the real line or scalings for use in
equilibrium conditions.

#### Fields

- `key::Symbol`: Parameter name. Should conform to the guidelines
  established in the DSGE Style Guide.
- `value::T`: The parameter's steady-state value.
- `description::AbstractString`: Short description of the parameter's economic significance.
- `tex_label::AbstractString`: String for printing parameter name to LaTeX.
"""
type SteadyStateParameter{T} <: AbstractParameter{T}
    key::Symbol
    value::T
    description::AbstractString
    tex_label::AbstractString
end

hasprior(p::Parameter) = !isnull(p.prior)

typealias NullableOrPrior Union{NullablePrior, ContinuousUnivariateDistribution}

# We want to use value field from UnscaledParameters and
# SteadyStateParameters in computation, so we alias their union here.
typealias UnscaledOrSteadyState Union{UnscaledParameter, SteadyStateParameter}

"""
```
ParamBoundsError <: Exception
```

A `ParamBoundsError` is thrown upon an attempt to assign a parameter value that is not
between `valuebounds`.
"""
type ParamBoundsError <: Exception
    msg::AbstractString
end
ParamBoundsError() = ParamBoundsError("Value not between valuebounds")
Base.showerror(io::IO, ex::ParamBoundsError) = print(io, ex.msg)

"""
```
parameter{T,U<:Transform}(key::Symbol, value::T, valuebounds = (value,value),
                          transform_parameterization = (value,value),
                          transform = Untransformed(), prior = NullablePrior(),
                          fixed = true, scaling::Function = identity, description = "",
                          tex_label::AbstractString = "")
```

By default, returns a fixed `UnscaledParameter` object with key `key`
and value `value`. If `scaling` is given, a `ScaledParameter` object
is returned.
"""

function parameter{T,U<:Transform}(key::Symbol,
                                   value::T,
                                   valuebounds::Interval{T} = (value,value),
                                   transform_parameterization::Interval{T} = (value,value),
                                   transform::U             = DSGE.Untransformed(),
                                   prior::NullableOrPrior   = NullablePrior();
                                   fixed::Bool              = true,
                                   scaling::Function        = identity,
                                   description::AbstractString = "",
                                   tex_label::AbstractString = "")

    # If fixed=true, force bounds to match and leave prior as null.  We need to define new
    # variable names here because of lexical scoping.

    valuebounds_new = valuebounds
    transform_parameterization_new = transform_parameterization
    transform_new = transform
    U_new = U
    prior_new = prior

    if fixed
        transform_parameterization_new = (value,value)  # value is transformed already
        transform_new = DSGE.Untransformed()            # fixed priors should stay untransformed
        U_new = Untransformed

        if isa(transform, Untransformed)
            valuebounds_new = (value,value)
        end
    else
        transform_parameterization_new = transform_parameterization
    end

    # ensure that we have a Nullable{Distribution}, if not construct one
    prior_new = !isa(prior_new,NullablePrior) ? NullablePrior(prior_new) : prior_new

    if scaling == identity
        return UnscaledParameter{T,U_new}(key, value, valuebounds_new,
                                          transform_parameterization_new, transform_new,
                                          prior_new, fixed, description, tex_label)
    else
        return ScaledParameter{T,U_new}(key, value, scaling(value), valuebounds_new,
                                        transform_parameterization_new, transform_new,
                                        prior_new, fixed, scaling, description, tex_label)
    end
end

"""
```
SteadyStateParameter{T<:Number}(key::Symbol, value::T;
                                description::AbstractString = "",
                                tex_label::AbstractString = "")
```

SteadyStateParameter constructor with optional `description` and `tex_label` arguments.
"""
function SteadyStateParameter{T<:Number}(key::Symbol,
                                       value::T;
                                       description::AbstractString = "",
                                       tex_label::AbstractString = "")

    return SteadyStateParameter(key, value, description, tex_label)
end

"""
```
parameter{T<:Number,U<:Transform}(p::UnscaledParameter{T,U}, newvalue::T)
```

Returns an UnscaledParameter with value field equal to `newvalue`. If `p` is a fixed
parameter, it is returned unchanged.
"""
function parameter{T<:Number,U<:Transform}(p::UnscaledParameter{T,U}, newvalue::T)
    p.fixed && return p    # if the parameter is fixed, don't change its value
    a,b = p.valuebounds
    if !(a <= newvalue <= b)
        throw(ParamBoundsError("New value of $(string(p.key)) ($(newvalue)) is out of bounds ($(p.valuebounds))"))
    end
    UnscaledParameter{T,U}(p.key, newvalue, p.valuebounds, p.transform_parameterization,
                           p.transform, p.prior, p.fixed, p.description, p.tex_label)
end

"""
```
parameter{T<:Number,U<:Transform}(p::ScaledParameter{T,U}, newvalue::T)
```

Returns a ScaledParameter with value field equal to `newvalue` and scaledvalue field equal
to `p.scaling(newvalue)`. If `p` is a fixed parameter, it is returned unchanged.
"""
function parameter{T<:Number,U<:Transform}(p::ScaledParameter{T,U}, newvalue::T)
    p.fixed && return p    # if the parameter is fixed, don't change its value
    a,b = p.valuebounds
    if !(a <= newvalue <= b)
        throw(ParamBoundsError("New value of $(string(p.key)) ($(newvalue)) is out of bounds ($(p.valuebounds))"))
    end
    ScaledParameter{T,U}(p.key, newvalue, p.scaling(newvalue), p.valuebounds,
                         p.transform_parameterization, p.transform, p.prior, p.fixed,
                         p.scaling, p.description, p.tex_label)
end

function Base.show{T,U}(io::IO, p::Parameter{T,U})
    @printf io "%s\n" typeof(p)
    @printf io "(:%s)\n%s\n"      p.key p.description
    @printf io "LaTeX label: %s\n"     p.tex_label
    @printf io "-----------------------------\n"
    #@printf io "real value:        %+6f\n" transform_to_real_line(p)
    @printf io "unscaled, untransformed value:        %+6f\n" p.value
    isa(p,ScaledParameter) && @printf "scaled, untransformed value:        %+6f\n" p.scaledvalue
    #!isa(U(),DSGE.Untransformed) && @printf io "transformed value: %+6f\n" p.value

    if hasprior(p)
        @printf io "prior distribution:\n\t%s\n" get(p.prior)
    else
        @printf io "prior distribution:\n\t%s\n" "no prior"
    end

    @printf io "transformation for csminwel:\n\t%s" U()
    @printf io "parameter is %s\n" p.fixed ? "fixed" : "not fixed"
end

function Base.show{T}(io::IO, p::SteadyStateParameter{T})
    @printf io "%s\n" typeof(p)
    @printf io "(:%s)\n%s\n"      p.key p.description
    @printf io "LaTeX label: %s\n"     p.tex_label
    @printf io "-----------------------------\n"
    @printf io "value:        %+6f\n" p.value
end

"""
```
transform_to_model_space{T<:Number, U<:Transform}(p::Parameter{T,U}, x::T)
```

Transforms `x` from the real line to lie between `p.valuebounds` without updating `p.value`.
The transformations are defined as follows, where (a,b) = p.transform_parameterization and c
a scalar (default=1):

- Untransformed: `x`
- SquareRoot:    `(a+b)/2 + (b-a)/2 * c * x/sqrt(1 + c^2 * x^2)`
- Exponential:   `a + exp(c*(x-b))`
"""
transform_to_model_space{T}(p::Parameter{T,Untransformed}, x::T) = x
function transform_to_model_space{T}(p::Parameter{T,SquareRoot}, x::T)
    (a,b), c = p.transform_parameterization, one(T)
    (a+b)/2 + (b-a)/2*c*x/sqrt(1 + c^2 * x^2)
end
function transform_to_model_space{T}(p::Parameter{T,Exponential}, x::T)
    (a,b),c = p.transform_parameterization,one(T)
    a + exp(c*(x-b))
end

transform_to_model_space{T}(pvec::ParameterVector{T}, values::Vector{T}) = map(transform_to_model_space, pvec, values)

"""
```
transform_to_real_line{T<:Number, U<:Transform}(p::Parameter{T,U}, x::T = p.value)
```

Transforms `p.value` from model space (between `p.valuebounds`) to the real line, without updating
`p.value`. The transformations are defined as follows,
where (a,b) = p.transform_parameterization, c a scalar (default=1), and x = p.value:

- Untransformed: x
- SquareRoot:   (1/c)*cx/sqrt(1 - cx^2), where cx =  2 * (x - (a+b)/2)/(b-a)
- Exponential:   a + exp(c*(x-b))
"""
transform_to_real_line{T}(p::Parameter{T,Untransformed}, x::T = p.value) = x
function transform_to_real_line{T}(p::Parameter{T,SquareRoot}, x::T = p.value)
    (a,b), c = p.transform_parameterization, one(T)
    cx = 2 * (x - (a+b)/2)/(b-a)
    (1/c)*cx/sqrt(1 - cx^2)
end
function transform_to_real_line{T}(p::Parameter{T,Exponential}, x::T = p.value)
    (a,b),c = p.transform_parameterization,one(T)
    b + (1/c) * log(x-a)
end

transform_to_real_line{T}(pvec::ParameterVector{T}, values::Vector{T}) = map(transform_to_real_line, pvec, values)
transform_to_real_line{T}(pvec::ParameterVector{T}) = map(transform_to_real_line, pvec)

# define operators to work on parameters

Base.convert{T<:Number}(::Type{T}, p::UnscaledParameter)  = convert(T,p.value)
Base.convert{T<:Number}(::Type{T}, p::ScaledParameter)    = convert(T,p.scaledvalue)
Base.convert{T<:Number}(::Type{T}, p::SteadyStateParameter)  = convert(T,p.value)

Base.promote_rule{T<:Number,U<:Number}(::Type{AbstractParameter{T}}, ::Type{U}) = promote_rule(T,U)

for op in (:(Base.(:+)),
           :(Base.(:-)),
           :(Base.(:*)),
           :(Base.(:/)),
           :(Base.(:^)))

    @eval ($op)(p::UnscaledOrSteadyState, q::UnscaledOrSteadyState) = ($op)(p.value, q.value)
    @eval ($op)(p::UnscaledOrSteadyState, x::Integer)            = ($op)(p.value, x)
    @eval ($op)(p::UnscaledOrSteadyState, x::Number)            = ($op)(p.value, x)
    @eval ($op)(x::Number, p::UnscaledOrSteadyState)            = ($op)(x, p.value)

    @eval ($op)(p::ScaledParameter, q::ScaledParameter) = ($op)(p.scaledvalue, q.scaledvalue)
    @eval ($op)(p::ScaledParameter, x::Integer)            = ($op)(p.scaledvalue, x)
    @eval ($op)(p::ScaledParameter, x::Number)            = ($op)(p.scaledvalue, x)
    @eval ($op)(x::Number, p::ScaledParameter)            = ($op)(x, p.scaledvalue)

    @eval ($op)(p::ScaledParameter, q::UnscaledOrSteadyState) = ($op)(p.scaledvalue, q.value)
    @eval ($op)(p::UnscaledOrSteadyState, q::ScaledParameter) = ($op)(p.value, q.scaledvalue)
end

for f in (:(Base.exp),
          :(Base.log),
          :(Base.(:-)),
          :(Base.(:<)),
          :(Base.(:>)),
          :(Base.(:<=)),
          :(Base.(:>=)))

    @eval ($f)(p::UnscaledOrSteadyState) = ($f)(p.value)
    @eval ($f)(p::ScaledParameter) = ($f)(p.scaledvalue)

    @eval ($f)(p::UnscaledOrSteadyState, x::Number) = ($f)(p.value, x)
    @eval ($f)(p::ScaledParameter, x::Number) = ($f)(p.scaledvalue, x)
end

"""
```
update!{T}(pvec::ParameterVector{T}, values::Vector{T})
```

Update all parameters in `pvec` that are not fixed with
`values`. Length of `values` must equal length of `pvec`.
"""
# this function is optimised for speed
function update!{T}(pvec::ParameterVector{T}, values::Vector{T})
    @assert length(values) == length(pvec) "Length of input vector (=$(length(values))) must match length of parameter vector (=$(length(pvec)))"
    map!(parameter, pvec, pvec, values)
end

"""
```
update{T}(pvec::ParameterVector{T}, values::Vector{T})
```

Returns a copy of `pvec` where non-fixed parameter values are updated
to `values`. `pvec` remains unchanged. Length of `values` must
equal length of `pvec`.
"""
# define the non-mutating version like this because we need the type stability of map!
update{T}(pvec::ParameterVector{T}, values::Vector{T}) = update!(copy(pvec), values)

Distributions.pdf(p::AbstractParameter) = exp(logpdf(p))
# we want the unscaled value for ScaledParameters
Distributions.logpdf{T,U}(p::Parameter{T,U}) = logpdf(get(p.prior),p.value)

# this function is optimised for speed
function Distributions.logpdf{T}(pvec::ParameterVector{T})
	x = zero(T)
	@inbounds for i = 1:length(pvec)
        if hasprior(pvec[i])
    		x += logpdf(pvec[i])
        end
	end
	x
end

# calculate logpdf at new values, without needing to allocate a temporary array with update
function Distributions.logpdf{T}(pvec::ParameterVector{T}, values::Vector{T})
    @assert length(values) == length(pvec) "Length of input vector (=$(length(values))) must match length of parameter vector (=$(length(pvec)))"

    x = zero(T)
    @inbounds for i = 1:length(pvec)
        if hasprior(pvec[i])
            x += logpdf(parameter(pvec[i], values[i]))
        end
    end
    x
end

Distributions.pdf{T}(pvec::ParameterVector{T}) = exp(logpdf(pvec))
Distributions.pdf{T}(pvec::ParameterVector{T}, values::Vector{T}) = exp(logpdf(pvec, values))