text:
println("including tssm_common.jl for type ", T)

function clone(psi::WaveFunction)
    wave_function(psi.m)
end

function eigen_function!(psi::WaveFunction1D, k)
    to_frequency_space!(psi)
    u=get_data(psi,true)
    u[:]=0.0
    u[k]=1
    psi
end

function eigen_function!(psi::WaveFunction2D, k, l)
    to_frequency_space!(psi)
    u=get_data(psi,true)
    u[:,:]=0.0
    u[k,l]=1
    psi
end

function eigen_function!(psi::WaveFunction3D, k, l, m)
    to_frequency_space!(psi)
    u=get_data(psi,true)
    u[:,:,:]=0.0
    u[k,l,m]=1
    psi
end

#NONSEPARATED_EIGENVALUES 0=cartesian, 1=polar/cylindrical, 2=spherical

for (METHOD, SUF, COMPLEX_METHOD, DIM, NONSEPARATED_EIGENVALUES) in (
                 (:Fourier1D, :_fourier_1d, true, 1, 0 ),
                 (:Fourier2D, :_fourier_2d, true, 2, 0 ),
                 (:Fourier3D, :_fourier_3d, true, 3, 0 ),
                 (:FourierReal1D, :_fourier_real_1d, false, 1, 0 ),
                 (:FourierReal2D, :_fourier_real_2d, false, 2, 0 ),
                 (:FourierReal3D, :_fourier_real_3d, false, 3, 0 ),
                 (:FourierBessel2D, :_fourier_bessel_2d, true, 2, 1 ),
                 (:FourierBesselReal2D, :_fourier_bessel_real_2d, false, 2, 1 ),
                 (:BesselRotSym1D, :_bessel_rotsym_1d, true, 1,0),
                 (:BesselRotSymReal1D, :_bessel_rotsym_real_1d, false, 1, 0),
                 (:Schroedinger1D, :_schroedinger_1d, true, 1, 0 ),
                 (:Schroedinger2D, :_schroedinger_2d, true, 2, 0 ),
                 (:Schroedinger3D, :_schroedinger_3d, true, 3, 0 ),
                 (:SchroedingerReal1D, :_schroedinger_real_1d, false, 1, 0 ),
                 (:SchroedingerReal2D, :_schroedinger_real_2d, false, 2, 0 ),
                 (:SchroedingerReal3D, :_schroedinger_real_3d, false, 3, 0 ),
                 (:SchroedingerRotating2D, :_schroedinger_rotating_2d, true, 2, 0 ),
                 (:SchroedingerRotating3D, :_schroedinger_rotating_3d, true, 3, 0 ),
                 (:SchroedingerRotatingReal2D, :_schroedinger_rotating_real_2d, false, 2, 0 ),
                 (:SchroedingerRotatingReal3D, :_schroedinger_rotating_real_3d, false, 3, 0 ),
                 (:SchroedingerHermite1D, :_schroedinger_hermite_1d, true, 1, 0 ),
                 (:SchroedingerHermite2D, :_schroedinger_hermite_2d, true, 2, 0 ),
                 (:SchroedingerHermite3D, :_schroedinger_hermite_3d, true, 3, 0 ),
                 (:SchroedingerHermiteReal1D, :_schroedinger_hermite_real_1d, false, 1, 0 ),
                 (:SchroedingerHermiteReal2D, :_schroedinger_hermite_real_2d, false, 2, 0 ),
                 (:SchroedingerHermiteReal3D, :_schroedinger_hermite_real_3d, false, 3, 0 ),
                 (:SchroedingerGeneralizedLaguerre2D, :_schroedinger_gen_laguerre_2d, true, 2, 0 ),
                 (:SchroedingerGeneralizedLaguerreHermite3D, :_schroedinger_gen_laguerre_hermite_3d, true, 3, 0 ),
                 (:SchroedingerGeneralizedLaguerreReal2D, :_schroedinger_gen_laguerre_real_2d, false, 2, 0 ),
                 (:SchroedingerGeneralizedLaguerreHermiteReal3D, :_schroedinger_gen_laguerre_hermite_real_3d, false, 3, 0 ),
                )
    println("    ", METHOD)
    if T == :Float128
        SUF = Symbol(SUF, "_wf128")
    end
    WF = Symbol(:Wf,METHOD)
    @eval begin
        # wave function constructor

        #function ($WF){($T)}( m::($METHOD){($T)} )
        function ($WF)( m::($METHOD){($T)} )
            wf = ($WF){($T)}( ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"new_wf",SUF))), Ptr{Nothing}, (Ptr{Nothing},), m.m) , m)
            finalizer(x -> ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"finalize_wf",SUF))), Nothing, (Ptr{Nothing},), x.p), wf )
            wf
        end
     end
     @eval begin

        wave_function(m::($METHOD){($T)}) = ($WF)(m)

        function set_propagate_time_together_with_A!(m::($METHOD){$T}, flag::Bool)
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_propagate_time_together_with_A",SUF))), Nothing,
               (Ptr{Nothing}, Int32), m.m, flag)
        end

        function get_propagate_time_together_with_A(m::($METHOD){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_propagate_time_together_with_A",SUF))), Int32,
               (Ptr{Nothing},), m.m) == 1
        end

        function is_real_space(psi::($WF){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"is_real_space_wf",SUF))), Int32,
               (Ptr{Nothing},), psi.p) == 1
        end

        function is_frequency_space(psi::($WF){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"is_real_space_wf",SUF))), Int32,
               (Ptr{Nothing},), psi.p) != 1
        end

        function to_real_space!(psi::($WF){$T})
             ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"to_real_space_wf",SUF))), Nothing,
                    (Ptr{Nothing},), psi.p)
        end

        function to_frequency_space!(psi::($WF){$T})
             ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"to_frequency_space_wf",SUF))), Nothing,
                    (Ptr{Nothing},), psi.p)
        end

        function save(psi::($WF){$T}, filename::AbstractString)
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"save_wf",SUF))), Nothing,
                 (Ptr{Nothing}, Cstring, Int32,), psi.p, filename, length(filename))
        end

        function load!(psi::($WF){$T}, filename::AbstractString)
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"load_wf",SUF))), Nothing,
                 (Ptr{Nothing}, Cstring, Int32,), psi.p, filename, length(filename))
        end

        function copy!(target::($WF){$T}, source::($WF){$T})
           if target.m ≠ source.m
               error("source and target must belong to the same method")
           end
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"copy_wf",SUF))), Nothing,
                 (Ptr{Nothing}, Ptr{Nothing}), target.p, source.p )
        end

        function norm(psi::($WF){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"norm_wf",SUF))), ($T),
                 (Ptr{Nothing}, ), psi.p )
        end

        function norm_in_frequency_space(psi::($WF){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"norm_in_frequency_space_wf",SUF))), ($T),
                 (Ptr{Nothing}, ), psi.p )
        end

        function distance(psi1::($WF){$T}, psi2::($WF){$T})
           if psi1.m ≠ psi2.m
               error("psi1 and psi2 must belong to the same method")
           end
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"distance_wf",SUF))), ($T),
                 (Ptr{Nothing}, Ptr{Nothing}), psi1.p, psi2.p )
        end

        function normalize!(psi::($WF){$T})
           ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"normalize_wf",SUF))), ($T),
                 (Ptr{Nothing}, ), psi.p )
        end

    end #eval

    if NONSEPARATED_EIGENVALUES==0 # cartesian
        if DIM==1
            @eval begin
                function get_eigenvalues(m::($METHOD){$T}, unsafe_access::Bool=false)
                   dim = zeros(Int32, 1)
                   evp = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_eigenvalues",SUF))), Ptr{$T},
                         (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 1 )
                   ev = unsafe_wrap(Array, evp, dim[1], own=false)
                   if unsafe_access
                       return ev
                   else
                       return copy(ev)
                   end
                end

            end
        elseif DIM==2
            @eval begin
               function get_eigenvalues(m::($METHOD){$T}, unsafe_access::Bool=false)
                   dim = zeros(Int32, 1)
                   evp1 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_eigenvalues",SUF))), Ptr{$T},
                         (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 1 )
                   ev1 = unsafe_wrap(Array, evp1, dim[1], own=false)
                   evp2 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_eigenvalues",SUF))), Ptr{$T},
                         (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 2 )
                   ev2 = unsafe_wrap(Array, evp2, dim[1], own=false)
                   if unsafe_access
                       return ev1, ev2
                   else
                       return copy(ev1), copy(ev2)
                   end
                end

            end
        elseif DIM==3
            @eval begin
               function get_eigenvalues(m::($METHOD){$T}, unsafe_access::Bool=false)
                   dim = zeros(Int32, 1)
                   evp1 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_eigenvalues",SUF))), Ptr{$T},
                         (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 1 )
                   ev1 = unsafe_wrap(Array, evp1, dim[1], own=false)
                   evp2 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_eigenvalues",SUF))), Ptr{$T},
                         (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 2 )
                   ev2 = unsafe_wrap(Array, evp2, dim[1], own=false)
                   evp3 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_eigenvalues",SUF))), Ptr{$T},
                         (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 3 )
                   ev3 = unsafe_wrap(Array, evp3, dim[1], own=false)
                   if unsafe_access
                       return ev1, ev2, ev3
                   else
                       return copy(ev1), copy(ev2), copy(ev3)
                   end
                end

            end
        end
    elseif NONSEPARATED_EIGENVALUES==1 # polar/cylindrical
        if DIM==2
            @eval begin
                function get_eigenvalues(m::($METHOD){$T}, unsafe_access::Bool=false)
                   dim = zeros(Int32, 2)
                   evp = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_eigenvalues",SUF))), Ptr{$T},
                         (Ptr{Nothing}, Ptr{Int32}), m.m, dim )
                   ev = unsafe_wrap(Array, evp, (dim[1], dim[2]), own=false)
                   if unsafe_access
                       return ev
                   else
                       return copy(ev)
                   end
                end

            end
        elseif DIM==3
            # tbd
        end
    elseif NONSEPARATED_EIGENVALUES==2 # spherical
        if DIM==3
            # tbd
        end
    end

    if DIM==1
        @eval begin

            function get_nodes(m::($METHOD){$T})
               dim = zeros(Int32, 1)
               np = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nodes",SUF))), Ptr{$T},
                     (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 1 )
               n = unsafe_wrap(Array, np, dim[1], own=false)
               copy(n)
            end

        end
    elseif DIM==2
        @eval begin

            function get_nodes(m::($METHOD){$T})
               dim = zeros(Int32, 1)
               np1 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nodes",SUF))), Ptr{$T},
                     (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 1 )
               n1 = unsafe_wrap(Array, np1, dim[1], own=false)
               np2 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nodes",SUF))), Ptr{$T},
                     (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 2 )
               n2 = unsafe_wrap(Array, np2, dim[1], own=false)
               copy(n1), copy(n2)
            end

        end
    elseif DIM==3
        @eval begin

            function get_nodes(m::($METHOD){$T})
               dim = zeros(Int32, 1)
               np1 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nodes",SUF))), Ptr{$T},
                     (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 1 )
               n1 = unsafe_wrap(Array, np1, dim[1], own=false)
               np2 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nodes",SUF))), Ptr{$T},
                     (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 2 )
               n2 = unsafe_wrap(Array, np2, dim[1], own=false)
               np3 = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_nodes",SUF))), Ptr{$T},
                     (Ptr{Nothing}, Ptr{Int32}, Int32), m.m, dim, 3 )
               n3 = unsafe_wrap(Array, np3, dim[1], own=false)
               copy(n1), copy(n2), copy(n3)
            end

          end
     end

    if COMPLEX_METHOD
        @eval begin

            function save(psi::($WF){$T}, filename::AbstractString,
                         dset_name_real::AbstractString, dset_name_imag::AbstractString;
                         append::Bool=false)
   	         ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"save1_wf",SUF))), Nothing,
                 	(Ptr{Nothing}, Cstring, Int32, Cstring, Int32, Cstring, Int32, Int32),
	                  psi.p, filename, length(filename),
        	          dset_name_real, length(dset_name_real),
                	  dset_name_imag, length(dset_name_imag), append)
            end

            function load!(psi::($WF){$T}, filename::AbstractString,
                         dset_name_real::AbstractString, dset_name_imag::AbstractString)
   	         ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"load1_wf",SUF))), Nothing,
                 	(Ptr{Nothing}, Cstring, Int32, Cstring, Int32, Cstring, Int32),
	                  psi.p, filename, length(filename),
        	          dset_name_real, length(dset_name_real),
                	  dset_name_imag, length(dset_name_imag))
            end

            function set_time!(psi::($WF){$T}, t::Number)
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_time_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Complex{$T},), psi.p, t)
            end

            function get_time(psi::($WF){$T})
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_time_wf",SUF))), ($T),
                        (Ptr{Nothing},), psi.p)
            end

            function propagate_time!(psi::($WF){$T}, dt::Number)
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"propagate_time_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Complex{$T},), psi.p, dt)
            end

            function propagate_A!(psi::($WF){$T}, dt::Number)
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"propagate_A_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Complex{$T},), psi.p, dt)
            end

            function propagate_A_derivative!(this::($WF){$T}, other::($WF){$T},
                                 dt::Number)
               if this.m ≠ other.m
                   error("this and other must belong to the same method")
               end
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"propagate_A_derivative_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Ptr{Nothing}, Complex{$T}),
                         this.p, other.p, dt)
            end

            function add_apply_A!(this::($WF){$T}, other::($WF){$T},
                                 coefficient::Number=1.0)
               if this.m ≠ other.m
                   error("this and other must belong to the same method")
               end
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"add_apply_A_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Ptr{Nothing}, Complex{$T}),
                         this.p, other.p, coefficient)
            end

            function add_phi_A!(this::($WF){$T}, other::($WF){$T}, dt::Number, n::Integer,
                                 coefficient::Number=1.0)
               if this.m ≠ other.m
                   error("this and other must belong to the same method")
               end
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"add_phi_A_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Ptr{Nothing}, Complex{$T}, Cint, Complex{$T}),
                         this.p, other.p, dt, n, coefficient)
            end

            function scale!(psi::($WF){$T}, factor::Number)
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"scale_wf",SUF))), Nothing,
                     (Ptr{Nothing}, Complex{$T} ), psi.p, factor )
            end

            function axpy!(this::($WF){$T}, other::($WF){$T},
                                 factor::Number)
               if this.m ≠ other.m
                   error("this and other must belong to the same method")
               end
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"axpy_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Ptr{Nothing}, Complex{$T}),
                         this.p, other.p, factor)
            end

            function inner_product(psi1::($WF){$T}, psi2::($WF){$T})
                if psi1.m ≠ psi2.m
                    error("psi1 and psi2 must belong to the same method")
                end
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"inner_product_wf",SUF))), (Complex{$T}),
                    (Ptr{Nothing}, Ptr{Nothing}), psi1.p, psi2.p )
            end

        end # eval

        if DIM==1
            @eval begin
                function set!(psi::($WF){$T}, x::Number)
                    u = get_data(psi, true)
                    u[:] .= x
                end

                function set!(psi::($WF){$T}, f::Function)
                   rt = Base.return_types(f, (($T),))
                   if length(rt)==1 && rt[1]==($T)
                       f_c = cfunction(f, ($T), (($T),))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"rset_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}), psi.p, f_c )
                   elseif length(rt)==1 && rt[1]==Complex{$T}
                       f_c = cfunction(f, Complex{$T}, (($T),))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}), psi.p, f_c )
                   else
                       error("wrong return type of function")
                   end
                end

                function set!(psi::($WF){$T}, f::Function, t::Real)
                   rt = Base.return_types(f, (($T),($T),))
                   if length(rt)==1 && rt[1]==($T)
                       f_c = cfunction(f, ($T), (($T), ($T),))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"rset_t_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}, ($T)), psi.p, f_c, t )
                   elseif length(rt)==1 && rt[1]==Complex{$T}
                       f_c = cfunction(f, Complex{$T}, (($T), ($T),))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_t_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}, ($T)), psi.p, f_c, t )
                   else
                       error("wrong return type of function")
                   end
                end

                function get_data(psi::($WF){$T}, unsafe_access::Bool=false)
                   dims = zeros(Int32, 1)
                   up = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_data_wf",SUF))), Ptr{Complex{$T}},
                         (Ptr{Nothing}, Ptr{Int32}), psi.p, dims )
                   data = unsafe_wrap(Array, up, dims[1], own=false)
                   if unsafe_access
                      return data
                   else
                      return copy(data)
                   end
                end

                function evaluate(psi::($WF){$T}, x::Real)
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"evaluate_wf",SUF))), Complex{$T},
                             (Ptr{Nothing}, ($T)), psi.p, x)
                end

            end # eval
        elseif DIM==2
            @eval begin

                function set!(psi::($WF){$T}, x::Number)
                    u = get_data(psi, true)
                    u[:,:] .= x
                end

                function set!(psi::($WF){$T}, f::Function)
                   rt = Base.return_types(f, (($T),($T),))
                   if length(rt)==1 && rt[1]==($T)
                       f_c = cfunction(f, ($T), (($T),($T)))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"rset_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}), psi.p, f_c )
                   elseif length(rt)==1 && rt[1]==Complex{$T}
                       f_c = cfunction(f, Complex{$T}, (($T),($T)))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}), psi.p, f_c )
                   else
                       error("wrong return type of function")
                   end
                end

                function set!(psi::($WF){$T}, f::Function, t::Real)
                   rt = Base.return_types(f, (($T),($T),($T),))
                   if length(rt)==1 && rt[1]==($T)
                       f_c = cfunction(f, ($T), (($T), ($T),($T)))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"rset_t_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}, ($T)), psi.p, f_c, t )
                   elseif length(rt)==1 && rt[1]==Complex{$T}
                       f_c = cfunction(f, Complex{$T}, (($T), ($T),($T)))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_t_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}, ($T)), psi.p, f_c, t )
                   else
                       error("wrong return type of function")
                   end
                end

                function get_data(psi::($WF){$T}, unsafe_access::Bool=false)
                   dims = zeros(Int32, 2)
                   up = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_data_wf",SUF))), Ptr{Complex{$T}},
                         (Ptr{Nothing}, Ptr{Int32}), psi.p, dims )
                   data = unsafe_wrap(Array, up, (dims[1], dims[2]), own=false)
                   if unsafe_access
                      return data
                   else
                      return copy(data)
                   end
                end

                function evaluate(psi::($WF){$T}, x::Real, y::Real)
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"evaluate_wf",SUF))), Complex{$T},
                             (Ptr{Nothing}, ($T),($T)), psi.p, x, y)
                end

            end # eval
        elseif DIM==3
            @eval begin
                function set!(psi::($WF){$T}, x::Number)
                    u = get_data(psi, true)
                    u[:,:,:] .= x
                end

                function set!(psi::($WF){$T}, f::Function)
                   rt = Base.return_types(f, (($T),($T),($T),))
                   if length(rt)==1 && rt[1]==($T)
                       f_c = cfunction(f, ($T), (($T),($T),($T)))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"rset_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}), psi.p, f_c )
                   elseif length(rt)==1 && rt[1]==Complex{$T}
                       f_c = cfunction(f, Complex{$T}, (($T),($T),($T)))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}), psi.p, f_c )
                   else
                       error("wrong return type of function")
                   end
                end

                function set!(psi::($WF){$T}, f::Function, t::Real)
                   rt = Base.return_types(f, (($T),($T),($T),($T)))
                   if length(rt)==1 && rt[1]==($T)
                       f_c = cfunction(f, ($T), (($T), ($T),($T),($T)))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"rset_t_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}, ($T)), psi.p, f_c, t )
                   elseif length(rt)==1 && rt[1]==Complex{$T}
                       f_c = cfunction(f, Complex{$T}, (($T), ($T),($T),($T)))
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_t_wf",SUF))), Nothing,
                             (Ptr{Nothing}, Ptr{Nothing}, ($T)), psi.p, f_c, t )
                   else
                       error("wrong return type of function")
                   end
                end

               function get_data(psi::($WF){$T}, unsafe_access::Bool=false)
                   dims = zeros(Int32, 3)
                   up = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_data_wf",SUF))), Ptr{Complex{$T}},
                         (Ptr{Nothing}, Ptr{Int32}), psi.p, dims )
                   data = unsafe_wrap(Array, up, (dims[1], dims[2], dim[3]), own=false)
                   if unsafe_access
                      return data
                   else
                      return copy(data)
                   end
                end

                function evaluate(psi::($WF){$T}, x::Real, y::Real, z::Real)
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"evaluate_wf",SUF))), Complex{$T},
                             (Ptr{Nothing}, ($T),($T),($T)), psi.p, x, y,z)
                end

            end # eval
        end

    else
        @eval begin

            function save(psi::($WF){$T}, filename::AbstractString,
                         dset_name::AbstractString;
			 append::Bool=false)
   	         ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"save1_wf",SUF))), Nothing,
                 	(Ptr{Nothing}, Cstring, Int32, Cstring, Int32),
	                  psi.p, filename, length(filename), dset_name, length(dset_name))
            end

            function load!(psi::($WF){$T}, filename::AbstractString,
                         dset_name::AbstractString)
   	         ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"load1_wf",SUF))), Nothing,
                 	(Ptr{Nothing}, Cstring, Int32, Cstring, Int32),
	                  psi.p, filename, length(filename), dset_name, length(dset_name))
            end

            function set_time!(psi::($WF){$T}, t::Real)
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_time_wf",SUF))), Nothing,
                        (Ptr{Nothing}, $T,), psi.p, t)
            end

            function get_time(psi::($WF){$T})
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_time_wf",SUF))), ($T),
                        (Ptr{Nothing},), psi.p)
            end

            function propagate_time!(psi::($WF){$T}, dt::Real)
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"propagate_time_wf",SUF))), Nothing,
                        (Ptr{Nothing}, $T,), psi.p, dt)
            end

            function propagate_A!(psi::($WF){$T}, dt::Real)
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"propagate_A_wf",SUF))), Nothing,
                        (Ptr{Nothing}, ($T),), psi.p, dt)
            end

            function propagate_A_derivative!(this::($WF){$T}, other::($WF){$T},
                                 dt::Real)
               if this.m ≠ other.m
                   error("this and other must belong to the same method")
               end
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"propagate_A_derivative_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Ptr{Nothing}, ($T)),
                         this.p, other.p, dt)
            end

            function add_apply_A!(this::($WF){$T}, other::($WF){$T},
                                 coefficient::Real=1.0)
               if this.m ≠ other.m
                   error("this and other must belong to the same method")
               end
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"add_apply_A_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Ptr{Nothing}, ($T)),
                         this.p, other.p, coefficient)
            end

            function add_phi_A!(this::($WF){$T}, other::($WF){$T}, dt::Real, n::Integer,
                                 coefficient::Real=1.0)
               if this.m ≠ other.m
                   error("this and other must belong to the same method")
               end
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"add_phi_A_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Ptr{Nothing}, ($T), Cint, ($T)),
                         this.p, other.p, dt, n, coefficient)
            end

            function scale!(psi::($WF){$T}, factor::Real)
               ccall(  Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"scale_wf",SUF))) , Nothing,
                     (Ptr{Nothing}, ($T) ), psi.p, factor )
            end

            function axpy!(this::($WF){$T}, other::($WF){$T},
                                 factor::Real)
               if this.m ≠ other.m
                   error("this and other must belong to the same method")
               end
               ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"axpy_wf",SUF))), Nothing,
                        (Ptr{Nothing}, Ptr{Nothing}, ($T)),
                         this.p, other.p, factor)
            end

            function inner_product(psi1::($WF){$T}, psi2::($WF){$T})
                if psi1.m ≠ psi2.m
                    error("psi1 and psi2 must belong to the same method")
                end
                ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"inner_product_wf",SUF))), ($T),
                    (Ptr{Nothing}, Ptr{Nothing}), psi1.p, psi2.p )
            end

        end # eval

        if DIM==1
            @eval begin

                function set!(psi::($WF){$T}, x::Real)
                    u = get_data(psi, true)
                    u[:] .= x
                end

                function set!(psi::($WF){$T}, f::Function)
                   f_c = cfunction_check_return_type(f, ($T), (($T),))
                   ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_wf",SUF))), Nothing,
                         (Ptr{Nothing}, Ptr{Nothing}), psi.p, f_c )
                end

                function set!(psi::($WF){$T}, f::Function, t::Real)
                   f_c = cfunction_check_return_type(f, ($T), (($T), ($T),))
                   ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_t_wf",SUF))), Nothing,
                         (Ptr{Nothing}, Ptr{Nothing}, ($T)), psi.p, f_c, t )
                end

                function get_data(psi::($WF){$T}, unsafe_access::Bool=false)
                   dims = zeros(Int32, 1)
                   up = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_data_wf",SUF))), Ptr{$T},
                         (Ptr{Nothing}, Ptr{Int32}), psi.p, dims )
                   data = unsafe_wrap(Array, up, dims[1], own=false)
                   if unsafe_access
                      return data
                   else
                      return copy(data)
                   end
                end

                function evaluate(psi::($WF){$T}, x::Real)
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"evaluate_wf",SUF))), ($T),
                             (Ptr{Nothing}, ($T)), psi.p, x)
                end

            end # eval
        elseif DIM==2
            @eval begin
                function set!(psi::($WF){$T}, x::Real)
                    u = get_data(psi, true)
                    u[:,:] .= x
                end

                function set!(psi::($WF){$T}, f::Function)
                   f_c = cfunction_check_return_type(f, ($T), (($T),($T)))
                   ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_wf",SUF))), Nothing,
                         (Ptr{Nothing}, Ptr{Nothing}), psi.p, f_c )
                end

                function set!(psi::($WF){$T}, f::Function, t::Real)
                   f_c = cfunction_check_return_type(f, ($T), (($T), ($T),($T)))
                   ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_t_wf",SUF))), Nothing,
                         (Ptr{Nothing}, Ptr{Nothing}, ($T)), psi.p, f_c, t )
                end

                function get_data(psi::($WF){$T}, unsafe_access::Bool=false)
                   dims = zeros(Int32, 2)
                   up = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_data_wf",SUF))), Ptr{$T},
                         (Ptr{Nothing}, Ptr{Int32}), psi.p, dims )
                   data = unsafe_wrap(Array, up, (dims[1], dims[2]), own=false)
                   if unsafe_access
                      return data
                   else
                      return copy(data)
                   end
                end

                function evaluate(psi::($WF){$T}, x::Real, y::Real)
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"evaluate_wf",SUF))), ($T),
                             (Ptr{Nothing}, ($T),($T)), psi.p, x, y)
                end

            end # eval
        elseif DIM==3
            @eval begin
                function set!(psi::($WF){$T}, x::Real)
                    u = get_data(psi, true)
                    u[:,:,:] .= x
                end

                function set!(psi::($WF){$T}, f::Function)
                   f_c = cfunction_check_return_type(f, ($T), (($T),($T),($T)))
                   ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_wf",SUF))), Nothing,
                         (Ptr{Nothing}, Ptr{Nothing}), psi.p, f_c )
                end

                function set!(psi::($WF){$T}, f::Function, t::Real)
                   f_c = cfunction_check_return_type(f, ($T), (($T), ($T),($T),($T)))
                   ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"set_t_wf",SUF))), Nothing,
                         (Ptr{Nothing}, Ptr{Nothing}, ($T)), psi.p, f_c, t )
                end

                function get_data(psi::($WF){$T}, unsafe_access::Bool=false)
                   dims = zeros(Int32, 3)
                   up = ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"get_data_wf",SUF))), Ptr{$T},
                         (Ptr{Nothing}, Ptr{Int32}), psi.p, dims )
                   data = unsafe_wrap(Array, up, (dims[1], dims[2], dim[3]), own=false)
                   if unsafe_access
                      return data
                   else
                      return copy(data)
                   end
                end

                function evaluate(psi::($WF){$T}, x::Real, y::Real, z::Real)
                       ccall( Libdl.dlsym(($TSSM_HANDLE), $(string(PRE,"evaluate_wf",SUF))), ($T),
                             (Ptr{Nothing}, ($T),($T),($T)), psi.p, x, y, z)
                end

            end # eval
        end
    end

end #for


text:
module TaskData
using Sessions, DataFrames, JLD
export stimulusframe, sessinfo, ISI

const ISI = [0.0, 0.2, 0.5, 0.8]

function stimulusframe(trials::DataFrame)
    trial_stimuli = trials[:stimuli].data::Vector{Vector{Int}}
    trial_onsets = trials[:onsets].data::Vector{Vector{Float64}}
    trial_offsets = trials[:offsets].data::Vector{Vector{Float64}}
    trial_correct::Vector{Bool} = trials[:correct].data
    trial_isi = trials[:isi].data::Vector{Float64}
    trial_probe1 = trials[:probe1].data::Vector{Int}
    trial_probe2 = trials[:probe2].data::Vector{Int}

    stimulus = zeros(Int, sum(length, trial_stimuli))
    next_stimulus = zeros(Int, sum(length, trial_stimuli))
    prev_stimulus = zeros(Int, sum(length, trial_stimuli))
    trial_number = zeros(Int, length(stimulus))
    pos = zeros(Int, length(stimulus))
    onset = fill(NaN, length(stimulus))
    offset = fill(NaN, length(stimulus))
    isi = fill(NaN, length(stimulus))
    probed = fill(false, length(stimulus))
    correct = DataArray(Bool, length(stimulus))
    n = 0
    for i = 1:length(trial_stimuli)
        rg = n+1:n+length(trial_stimuli[i])
        trial_number[rg] = i
        pos[rg] = 1:length(rg)
        stimulus[rg] = trial_stimuli[i]
        prev_stimulus[rg[2:end]] = trial_stimuli[i][1:end-1]
        next_stimulus[rg[1:end-1]] = trial_stimuli[i][2:end]
        onset[rg] = trial_onsets[i]
        offset[rg] = trial_offsets[i]
        isi[rg] = trial_isi[i]
        pb = probed[rg] = (trial_stimuli[i] .== trial_probe1[i]) | (trial_stimuli[i] .== trial_probe2[i])
        @assert sum(pb) == 1
        correct[rg[findfirst(pb)]] = trial_correct[i]
        n = rg[end]
    end
    @assert n == length(stimulus)

    DataFrame(trial_number=trial_number, pos=pool(pos), stimulus=pool(stimulus),
              prev_stimulus=pool(prev_stimulus), next_stimulus=pool(next_stimulus),
              onset=onset, offset=offset, isi=pool(isi), probed=probed, correct=correct)
end

function sessinfo()
    patients = Int[]
    correct = Float64[]
    trials = Int[]
    stimuli = Int[]
    npertrial = Int[]
    maintenancetime = Tuple{Float64,Float64,Float64}[]
    snames = sessionnames()
    for sname in snames
        s = readsession(sname)
        push!(patients, parse(Int, sname[1:3]))
        push!(correct, mean(s.trials[:correct]))
        push!(trials, size(s.trials, 1))
        push!(stimuli, length(s.stimuli))
        push!(npertrial, length(s.trials[:stimuli][1]))
        maintenancetimes = s.trials[:probes] .- s.trials[:maint]
        push!(maintenancetime, (median(maintenancetimes), extrema(maintenancetimes)...))
    end
    DataFrame(session=snames, correct=correct, trials=trials, stimuli=stimuli,
              npertrial=npertrial, maintenancetime=maintenancetime)
end
end

ext: jl

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
@testset "AlgAss" begin
  include("AlgAss/AlgAss.jl")
  include("AlgAss/AlgGrp.jl")
  include("AlgAss/Elem.jl")
end


text:
# functions for computing the jacobians of the boundary conditions calculated
# in bc.jl

"""
  Computes the jacobian of the boundary condition terms computed by
  [`getBCFluxes`](@ref) and assembles them into the Jacobian

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * assembler: AssembleElementData
"""
function getBCFluxes_diff(mesh::AbstractMesh, sbp::AbstractOperator, eqn::EulerData,
                          opts, assembler::AssembleElementData)
  #get all the fluxes for all the boundary conditions and save them in eqn.bndryflux

  #println("mesh.bndry_funcs = ", mesh.bndry_funcs)
  for i=1:mesh.numBC
  #  println("computing flux for boundary condition ", i)
    functor_i = mesh.bndry_funcs[i]
    start_index = mesh.bndry_offsets[i]
    end_index = mesh.bndry_offsets[i+1]
    idx_range = start_index:(end_index - 1)
    bndry_facenums_i = sview(mesh.bndryfaces, start_index:(end_index - 1))

    calcBoundaryFlux_nopre_diff(mesh, sbp, eqn, functor_i, idx_range, bndry_facenums_i, assembler)
  end

  return nothing
end

"""
### EulerEquationMod.getBCFluxes_revm

Reverse mode of getBCFluxes w.r.t mesh metrics

**Arguments**

* mesh  : a mesh object
* sbp   : SBP operator object
* eqn   : an EulerData object
* opts  : options dictionary

"""
function getBCFluxes_revm(mesh::AbstractMesh, sbp::AbstractOperator, eqn::EulerData, opts)
  #get all the fluxes for all the boundary conditions and save them in eqn.bndryflux

#  fill!(mesh.dxidx_bndry_bar, 0.0)
  for i=1:mesh.numBC
    functor_i = mesh.bndry_funcs_revm[i]
    start_index = mesh.bndry_offsets[i]
    end_index = mesh.bndry_offsets[i+1] - 1
    idx_range = start_index:(end_index)
    bndry_facenums_i = sview(mesh.bndryfaces, start_index:end_index)

    # call the function that calculates the flux for this boundary condition
    # passing the functor into another function avoid type instability
    calcBoundaryFlux_nopre_revm(mesh, sbp, eqn, functor_i, idx_range, bndry_facenums_i)
  end

  writeBoundary(mesh, sbp, eqn, opts)

  return nothing
end

"""
  Reverse mode of [`getBCFluxes`](@ref) wrt the solution.  `eqn.res_bar` is the
  input and `eqn.q_bar` is updated with the result
"""
function getBCFluxes_revq(mesh::AbstractMesh, sbp::AbstractOperator, eqn::EulerData, opts)
  #get all the fluxes for all the boundary conditions and save them in eqn.bndryflux

#  fill!(mesh.dxidx_bndry_bar, 0.0)
  for i=1:mesh.numBC
    functor_i = mesh.bndry_funcs_revq[i]
    start_index = mesh.bndry_offsets[i]
    end_index = mesh.bndry_offsets[i+1] - 1
    idx_range = start_index:(end_index)
    bndry_facenums_i = sview(mesh.bndryfaces, start_index:end_index)

    # call the function that calculates the flux for this boundary condition
    # passing the functor into another function avoid type instability
    calcBoundaryFlux_nopre_revq(mesh, sbp, eqn, functor_i, idx_range, bndry_facenums_i)
  end

  writeBoundary(mesh, sbp, eqn, opts)

  return nothing
end

"""
  Computes the jacobian of the boundary condition terms for a given
  BC and assembles it into the Jacobian.  See [`calcBoundaryFlux_nopre`](@ref).
  Note that this function does not depend on any quatiites being precomputed.

  **Inputs**

   * mesh
   * sbp
   * eqn
   * functor: this is the regular BCType functor, not BCType_diff
   * idx_range: range of indices in mesh.bndryfaces that bndry_facenums
                came from
   * bndry_facenums: array of Boundary objects to compute the BC for
   * assembler: AssembleElementData needed to assemble the element level
                jacobian into the system jacobian

  **Implementation Notes**

  Currently this computes the jacobian of the flux with respect to the
  face via complex step (using eqn.params_complex).  This should be replaced
  with dual numbers eventually
"""
function calcBoundaryFlux_nopre_diff(mesh::AbstractDGMesh{Tmsh},
                          sbp::AbstractOperator, eqn::EulerData{Tsol1, Tres1},
                          functor::BCType, idx_range::UnitRange,
                          bndry_facenums::AbstractArray{Boundary,1},
                          assembler::AssembleElementData) where {Tsol1, Tres1, Tmsh}
  # calculate the jacobian of the boundary flux for the boundary condition
  # evaluated by the functor
  # because we are complex stepping the boundary condition and Tsol is usually
  # a real number, the Tsol and Tres of the eqn object are not the right
  # Tsol and Tres for this function

  nfaces = length(bndry_facenums)
  params = eqn.params_complex
  q_face = params.calc_face_integrals_data.q_faceL
  Tsol = eltype(q_face)
  Tres = promote_type(Tsol, Tmsh)

#  flux_face = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  aux_vars = Array{Tres}(1)
  flux_k = zeros(Tres, mesh.numDofPerNode)
  flux_jac = params.calc_face_integrals_data.flux_dotL
  res_jac = params.calc_face_integrals_data.res_jacLL

  h = 1e-20
  pert = Tsol(0, h)
  for i=1:nfaces  # loop over faces with this BC
    bndry_i = bndry_facenums[i]
    global_facenum = idx_range[i]

    # interpolate to face
    q_volume = sview(eqn.q, :, :, bndry_i.element)
    boundaryFaceInterpolate!(mesh.sbpface, bndry_i.face, q_volume, q_face)

    for j = 1:mesh.numNodesPerFace

      # get components
      # q = ro_sview(eqn.q_bndry, :, j, global_facenum)
      q_j = sview(q_face, :, j)
      coords = ro_sview(mesh.coords_bndry, :, j, global_facenum)
      nrm_xy = ro_sview(mesh.nrm_bndry, :, j, global_facenum)

      bndry_node = BoundaryNode(bndry_i, i, j)
      # compute the jacobian of the flux wrt q_face
      for k=1:mesh.numDofPerNode
        q_j[k] += pert
        aux_vars[1] = calcPressure(params, q_j)

        functor(params, q_j, aux_vars, coords, nrm_xy, flux_k, bndry_node)

        for p=1:mesh.numDofPerNode
          flux_jac[p, k, j] = imag(flux_k[p])/h
        end

        q_j[k] -= pert
      end  # end loop k
    end  # end looop j

    fill!(res_jac, 0.0)
    boundaryFaceCombined_jac!(mesh.sbpface, bndry_i.face, flux_jac, res_jac,
                               SummationByParts.Subtract())

    if eqn.params.use_Minv == 1
      # multiply by Minv if needed
      for q=1:mesh.numNodesPerElement
        for p=1:mesh.numNodesPerElement
          val = mesh.jac[p, bndry_i.element]/sbp.w[p]  # entry in Minv
          @simd for m=1:mesh.numDofPerNode
            @simd for n=1:mesh.numDofPerNode
              res_jac[n, m, p, q] *= val
            end
          end
        end
      end
    end

    assembleBoundary(assembler, mesh.sbpface, mesh, bndry_i, res_jac)

  end  # end loop i

  return nothing
end

"""
  Reverse mode wrt metrics of [`calcBoundaryFlux_nopre`](@ref)
"""
function calcBoundaryFlux_nopre_revm( mesh::AbstractDGMesh{Tmsh},
                          sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                          functor_bar::BCType_revm, idx_range::UnitRange,
                          bndry_facenums::AbstractArray{Boundary,1}) where {Tmsh,  Tsol, Tres}
  # calculate the boundary flux for the boundary condition evaluated by the
  # functor

  nfaces = length(bndry_facenums)
  q2 = zeros(Tsol, mesh.numDofPerNode)
  params = eqn.params
  flux_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  for i=1:nfaces  # loop over faces with this BC
    bndry_i = bndry_facenums[i]
    global_facenum = idx_range[i]

    res_bar_i = sview(eqn.res_bar, :, :, bndry_i.element)
    fill!(flux_face_bar, 0)
    boundaryFaceIntegrate_rev!(mesh.sbpface, bndry_i.face, flux_face_bar, res_bar_i,
                           SummationByParts.Subtract())

    for j = 1:mesh.numNodesPerFace

      # get components
      #TODO: this doesn't work if precompute_q_bndr == false ?
      q = ro_sview(eqn.q_bndry, :, j, global_facenum)
      # convert to conservative variables if needed
      convertToConservative(eqn.params, q, q2)
      aux_vars = ro_sview(eqn.aux_vars_bndry, :, j, global_facenum)
      coords = ro_sview(mesh.coords_bndry, :, j, global_facenum)
      coords_bar = sview(mesh.coords_bndry_bar, :, j, global_facenum)
      nrm_xy = ro_sview(mesh.nrm_bndry, :, j, global_facenum)
      nrm_bar = sview(mesh.nrm_bndry_bar, :, j, global_facenum)
      bndryflux_bar_i = sview(flux_face_bar, :, j)

      bndry_node = BoundaryNode(bndry_i, i, j)

      functor_bar(params, q2, aux_vars, coords, coords_bar, nrm_xy, nrm_bar,
                  bndryflux_bar_i, bndry_node)
    end

  end

  return nothing
end

"""
  Reverse mode wrt metrics of [`calcBoundaryFlux_nopre`](@ref)
"""
function calcBoundaryFlux_nopre_revq( mesh::AbstractDGMesh{Tmsh},
                          sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                          functor_bar::BCType_revq, idx_range::UnitRange,
                          bndry_facenums::AbstractArray{Boundary,1}) where {Tmsh,  Tsol, Tres}
  # calculate the boundary flux for the boundary condition evaluated by the
  # functor

  nfaces = length(bndry_facenums)
#  q2 = zeros(Tsol, mesh.numDofPerNode)
  q_face = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  params = eqn.params
  flux_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  for i=1:nfaces  # loop over faces with this BC
    bndry_i = bndry_facenums[i]
    global_facenum = idx_range[i]

    q_vol = sview(eqn.q, :, :, bndry_i.element)
    q_vol_bar = sview(eqn.q_bar, :, :, bndry_i.element)
    boundaryFaceInterpolate!(mesh.sbpface, bndry_i.face, q_vol, q_face)

    res_bar_i = sview(eqn.res_bar, :, :, bndry_i.element)
    fill!(flux_face_bar, 0)
    boundaryFaceIntegrate_rev!(mesh.sbpface, bndry_i.face, flux_face_bar,
                               res_bar_i, SummationByParts.Subtract())

    fill!(q_face_bar, 0)

    for j = 1:mesh.numNodesPerFace
      # get components
      q_j = sview(q_face, :, j)
      q_bar_j = sview(q_face_bar, :, j)
      # convert to conservative variables if needed
      #convertToConservative(eqn.params, q, q2)
      aux_vars = ro_sview(eqn.aux_vars_bndry, :, j, global_facenum)
      coords = ro_sview(mesh.coords_bndry, :, j, global_facenum)
      nrm_xy = ro_sview(mesh.nrm_bndry, :, j, global_facenum)
      bndryflux_bar_i = sview(flux_face_bar, :, j)

      bndry_node = BoundaryNode(bndry_i, i, j)

      functor_bar(params, q_j, q_bar_j, aux_vars, coords, nrm_xy,
                  bndryflux_bar_i, bndry_node)
    end

    boundaryFaceInterpolate_rev!(mesh.sbpface, bndry_i.face, q_vol_bar, q_face_bar)
  end

  return nothing
end


text:
using CVChannel
using Test

"""
This script verifies the multiplicativity of the tensor product between
the Siddhu channel and the generalized amplitude damping channel (GADC).
The GADC is defined in https://arxiv.org/abs/2107.13486
"""

print("\nHere we verify that the Siddhu channel is multiplicative with GADC.")
@testset "Verifying Multiplicativity of Siddhu Channel with GDAC" begin
    sid_scan_range =[0:0.1:0.5;]
    gad_scan_range = [0:0.1:1;]
    sid_ctr, gad_ctr_1, gad_ctr_2 = 1,1,1
    max_val = 0
    is_mult = true
    for s_id in sid_scan_range
        println("---Now scanning for s =", s_id, ".-----")
        for p_id in gad_scan_range
            println("Now scanning for p =", p_id,".")
            target_val = 2*(1+sqrt(1-p_id))
            for n_id in gad_scan_range
                symgadchan(X) = GADChannel(X,p_id,n_id)
                sidchan(X) = siddhuChannel(X,s_id)

                gad_chan = Choi(symgadchan,2,2)
                sid_chan= Choi(sidchan,3,3)

                par_cv, = pptCV(parChoi(sid_chan, gad_chan), :dual)

                non_mult = par_cv - target_val
                !isapprox(non_mult,0,atol=5e-6) ? is_mult = false : nothing
                max_val < non_mult ? max_val = non_mult : nothing

                gad_ctr_2 += 1
            end
            gad_ctr_1 += 1
            gad_ctr_2 = 1
        end
        sid_ctr += 1
        gad_ctr_1 = 1
    end
    @test is_mult
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
"""
# Usage

    solve(prob::SingleTermFODEProblem, h, GL())
Grunwald Letnikov method for fractional ordinary differential equations

```tex
@INPROCEEDINGS{8742063,
author={Clemente-López, D. and Muñoz-Pacheco, J. M. and Félix-Beltrán, O. G. and Volos, C.},
booktitle={2019 8th International Conference on Modern Circuits and Systems Technologies (MOCAST)},
title={Efficient Computation of the Grünwald-Letnikov Method for ARM-Based Implementations of Fractional-Order Chaotic Systems},
year={2019},
doi={10.1109/MOCAST.2019.8742063}}
```
"""
struct GL <: FractionalDiffEqAlgorithm end

function solve(FODE::SingleTermFODEProblem, h, ::GL)
    @unpack f, α, u0, tspan = FODE
    t0 = tspan[1]; T = tspan[2]
    N::Int = floor(Int, (T-t0)/h)+1
    c = zeros(Float64, N)

    cp::Float64 = 1.0
    for j = 1:N
        c[j] = (1-(1+α)/j)*cp
        cp = c[j]
    end

    # Initialization
    y = zeros(Float64, N)
    y[1] = u0

    @fastmath @inbounds @simd for i = 2:N
        right = 0
        @fastmath @inbounds @simd for j=1:i-1
            right += c[j]*y[i-j]
        end
        y[i] = f(t0+(i-1)*h, y[i-1])*h^α - right
    end
    return FODESolution(collect(t0:h:T), y)
end


text:
#!/usr/local/bin/julia

using ArgParse

function parse_commandline()
    settings = ArgParseSettings()

    @add_arg_table settings begin
        "--alpha"
            help = "begin with alpha"
        "--omega", "-o"
	    help = "end with omega"
            arg_type = Bool
            default = true
	"--end"
	    help = "end the universe"
	    action = :store_true
        "arg"
            help = "positional arguments"
            required = true
    end

    return parse_args(settings)
end

"""
./args.jl --alpha A -o true goof

"""
function main()
    parsed_args = parse_commandline()
    println("args:")
    for (arg,val) in parsed_args
	    println("  [$arg]  =>  [$val]")
    end
end

main()

max_stars_repo_name: rabarar/juliaSnippets

max_issues_repo_name: rabarar/juliaSnippets

max_forks_repo_name: rabarar/juliaSnippets


text:
# estimate GARCH(1,1) with AR(1) in mean
using Statistics, DelimitedFiles, Econometrics
function main()
include("garch11.jl")
# weekly close price of NSYE, data provided with GRETL
data = readdlm("nysewk.txt")
# compute weekly percentage growth
y = 100.0 * log.(data[2:end] ./ data[1:end-1])
# Constrained maximization of logL
# note that the objective has a minus sign in front, as fmincon
# minimizes, but we want to maximize the logL
# get constrained estimates to use as input for MLE
thetastart = [mean(y); 0.0; var(y); 0.1; 0.1]
obj = theta -> -sum(garch11(theta, y))
thetahat, logL, junk  = fmincon(obj, thetastart, [], [], [-Inf, -1.0, 0.0, 0.0, 0.0], [Inf, 1.0, Inf, 1.0, 1.0])
# now do MLE to see formatted results.
# NOTE TO SELF: this won't work when the constraints are binding
# should add a method to mle for constrained problems.
model = theta -> garch11(theta, y)
thetahat, logL, junk, converged = mleresults(model, thetahat, "GARCH(1,1) example")
return
end
main()

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
# # Ginzburg-Landau model energy minimization

# ![landau_orig.png](landau_orig.png)

# Original

# ![landau_opt.png](landau_opt.png)

# Optimized

# In this example a basic Ginzburg-Landau model is solved.
# This example gives an idea of how the API together with ForwardDiff can be leveraged to
# performantly solve non standard problems on a FEM grid.
# A large portion of the code is there only for performance reasons,
# but since this usually really matters and is what takes the most time to optimize,
# it is included.

# The key to using a method like this for minimizing a free energy function directly,
# rather than the weak form, as is usually done with FEM, is to split up the
# gradient and Hessian calculations.
# This means that they are performed for each cell seperately instead of for the
# grid as a whole.

using ForwardDiff
import ForwardDiff: GradientConfig, HessianConfig, Chunk
using Ferrite
using Optim, LineSearches
using SparseArrays
using Tensors
using Base.Threads
# ## Energy terms
# ### 4th order Landau free energy
function Fl(P::Vec{3, T}, α::Vec{3}) where T
    P2 = Vec{3, T}((P[1]^2, P[2]^2, P[3]^2))
    return (α[1] * sum(P2) +
           α[2] * (P[1]^4 + P[2]^4 + P[3]^4)) +
           α[3] * ((P2[1] * P2[2]  + P2[2]*P2[3]) + P2[1]*P2[3])
end
# ### Ginzburg free energy
@inline Fg(∇P, G) = 0.5(∇P ⊡ G) ⊡ ∇P
# ### GL free energy
F(P, ∇P, params)  = Fl(P, params.α) + Fg(∇P, params.G)

# ### Parameters that characterize the model
struct ModelParams{V, T}
    α::V
    G::T
end

# ### ThreadCache
# This holds the values that each thread will use during the assembly.
struct ThreadCache{CV, T, DIM, F <: Function, GC <: GradientConfig, HC <: HessianConfig}
    cvP              ::CV
    element_indices  ::Vector{Int}
    element_dofs     ::Vector{T}
    element_gradient ::Vector{T}
    element_hessian  ::Matrix{T}
    element_coords   ::Vector{Vec{DIM, T}}
    element_potential::F
    gradconf         ::GC
    hessconf         ::HC
end
function ThreadCache(dpc::Int, nodespercell, cvP::CellValues{DIM, T}, modelparams, elpotential) where {DIM, T}
    element_indices  = zeros(Int, dpc)
    element_dofs     = zeros(dpc)
    element_gradient = zeros(dpc)
    element_hessian  = zeros(dpc, dpc)
    element_coords   = zeros(Vec{DIM, T}, nodespercell)
    potfunc          = x -> elpotential(x, cvP, modelparams)
    gradconf         = GradientConfig(potfunc, zeros(dpc), Chunk{12}())
    hessconf         = HessianConfig(potfunc, zeros(dpc), Chunk{12}())
    return ThreadCache(cvP, element_indices, element_dofs, element_gradient, element_hessian, element_coords, potfunc, gradconf, hessconf)
end

# ## The Model
# everything is combined into a model.
mutable struct LandauModel{T, DH <: DofHandler, CH <: ConstraintHandler, TC <: ThreadCache}
    dofs          ::Vector{T}
    dofhandler    ::DH
    boundaryconds ::CH
    threadindices ::Vector{Vector{Int}}
    threadcaches  ::Vector{TC}
end

function LandauModel(α, G, gridsize, left::Vec{DIM, T}, right::Vec{DIM, T}, elpotential) where {DIM, T}
    grid = generate_grid(Tetrahedron, gridsize, left, right)
    threadindices = Ferrite.create_coloring(grid)

    qr  = QuadratureRule{DIM, RefTetrahedron}(2)
    cvP = CellVectorValues(qr, Lagrange{DIM, RefTetrahedron, 1}())

    dofhandler = DofHandler(grid)
    push!(dofhandler, :P, 3)
    close!(dofhandler)

    dofvector = zeros(ndofs(dofhandler))
    startingconditions!(dofvector, dofhandler)
    boundaryconds = ConstraintHandler(dofhandler)
    #boundary conditions can be added but aren't necessary for optimization
    #add!(boundaryconds, Dirichlet(:P, getfaceset(grid, "left"), (x, t) -> [0.0,0.0,0.53], [1,2,3]))
    #add!(boundaryconds, Dirichlet(:P, getfaceset(grid, "right"), (x, t) -> [0.0,0.0,-0.53], [1,2,3]))
    close!(boundaryconds)
    update!(boundaryconds, 0.0)

    apply!(dofvector, boundaryconds)

    hessian = create_sparsity_pattern(dofhandler)
    dpc = ndofs_per_cell(dofhandler)
    cpc = length(grid.cells[1].nodes)
    caches = [ThreadCache(dpc, cpc, copy(cvP), ModelParams(α, G), elpotential) for t=1:nthreads()]
    return LandauModel(dofvector, dofhandler, boundaryconds, threadindices, caches)
end

# utility to quickly save a model
function Ferrite.vtk_save(path, model, dofs=model.dofs)
    vtkfile = vtk_grid(path, model.dofhandler)
    vtk_point_data(vtkfile, model.dofhandler, dofs)
    vtk_save(vtkfile)
end

# ## Assembly
# This macro defines most of the assembly step, since the structure is the same for
# the energy, gradient and Hessian calculations.
macro assemble!(innerbody)
    esc(quote
        dofhandler = model.dofhandler
        for indices in model.threadindices
            @threads for i in indices
                cache     = model.threadcaches[threadid()]
                eldofs    = cache.element_dofs
                nodeids   = dofhandler.grid.cells[i].nodes
                for j=1:length(cache.element_coords)
                    cache.element_coords[j] = dofhandler.grid.nodes[nodeids[j]].x
                end
                reinit!(cache.cvP, cache.element_coords)

                celldofs!(cache.element_indices, dofhandler, i)
                for j=1:length(cache.element_dofs)
                    eldofs[j] = dofvector[cache.element_indices[j]]
                end
                $innerbody
            end
        end
    end)
end

# This calculates the total energy calculation of the grid
function F(dofvector::Vector{T}, model) where T
    outs = fill(zero(T), nthreads())
    @assemble! begin
        outs[threadid()] += cache.element_potential(eldofs)
    end
    return sum(outs)
end

# The gradient calculation for each dof
function ∇F!(∇f::Vector{T}, dofvector::Vector{T}, model::LandauModel{T}) where T
    fill!(∇f, zero(T))
    @assemble! begin
        ForwardDiff.gradient!(cache.element_gradient, cache.element_potential, eldofs, cache.gradconf)
        @inbounds assemble!(∇f, cache.element_indices, cache.element_gradient)
    end
end

# The Hessian calculation for the whole grid
function ∇²F!(∇²f::SparseMatrixCSC, dofvector::Vector{T}, model::LandauModel{T}) where T
    assemblers = [start_assemble(∇²f) for t=1:nthreads()]
    @assemble! begin
        ForwardDiff.hessian!(cache.element_hessian, cache.element_potential, eldofs, cache.hessconf)
        @inbounds assemble!(assemblers[threadid()], cache.element_indices, cache.element_hessian)
    end
end

# We can also calculate all things in one go!
function calcall(∇²f::SparseMatrixCSC, ∇f::Vector{T}, dofvector::Vector{T}, model::LandauModel{T}) where T
    outs = fill(zero(T), nthreads())
    fill!(∇f, zero(T))
    assemblers = [start_assemble(∇²f, ∇f) for t=1:nthreads()]
    @assemble! begin
        outs[threadid()] += cache.element_potential(eldofs)
        ForwardDiff.hessian!(cache.element_hessian, cache.element_potential, eldofs, cache.hessconf)
        ForwardDiff.gradient!(cache.element_gradient, cache.element_potential, eldofs, cache.gradconf)
        @inbounds assemble!(assemblers[threadid()], cache.element_indices, cache.element_gradient, cache.element_hessian)
    end
    return sum(outs)
end

# ## Minimization
# Now everything can be combined to minimize the energy, and find the equilibrium
# configuration.
function minimize!(model; kwargs...)
    dh = model.dofhandler
    dofs = model.dofs
    ∇f = fill(0.0, length(dofs))
    ∇²f = create_sparsity_pattern(dh)
    function g!(storage, x)
        ∇F!(storage, x, model)
        apply_zero!(storage, model.boundaryconds)
    end
    function h!(storage, x)
        ∇²F!(storage, x, model)
        #apply!(storage, model.boundaryconds)
    end
    f(x) = F(x, model)

    od = TwiceDifferentiable(f, g!, h!, model.dofs, 0.0, ∇f, ∇²f)

    # this way of minimizing is only beneficial when the initial guess is completely off,
    # then a quick couple of ConjuageGradient steps brings us easily closer to the minimum.
    # res = optimize(od, model.dofs, ConjugateGradient(linesearch=BackTracking()), Optim.Options(show_trace=true, show_every=1, g_tol=1e-20, iterations=10))
    # model.dofs .= res.minimizer
    # to get the final convergence, Newton's method is more ideal since the energy landscape should be almost parabolic
    #+
    res = optimize(od, model.dofs, Newton(linesearch=BackTracking()), Optim.Options(show_trace=true, show_every=1, g_tol=1e-20))
    model.dofs .= res.minimizer
    return res
end

# ## Testing it
# This calculates the contribution of each element to the total energy,
# it is also the function that will be put through ForwardDiff for the gradient and Hessian.
function element_potential(eldofs::AbstractVector{T}, cvP, params) where T
    energy = zero(T)
    for qp=1:getnquadpoints(cvP)
        P  = function_value(cvP, qp, eldofs)
        ∇P = function_gradient(cvP, qp, eldofs)
        energy += F(P, ∇P, params) * getdetJdV(cvP, qp)
    end
    return energy
end

# now we define some starting conditions
function startingconditions!(dofvector, dofhandler)
    for cell in CellIterator(dofhandler)
        globaldofs = celldofs(cell)
        it = 1
        for i=1:3:length(globaldofs)
            dofvector[globaldofs[i]]   = -2.0
            dofvector[globaldofs[i+1]] = 2.0
            dofvector[globaldofs[i+2]] = -2.0tanh(cell.coords[it][1]/20)
            it += 1
        end
    end
end

δ(i, j) = i == j ? one(i) : zero(i)
V2T(p11, p12, p44) = Tensor{4, 3}((i,j,k,l) -> p11 * δ(i,j)*δ(k,l)*δ(i,k) + p12*δ(i,j)*δ(k,l)*(1 - δ(i,k)) + p44*δ(i,k)*δ(j,l)*(1 - δ(i,j)))

G = V2T(1.0e2, 0.0, 1.0e2)
α = Vec{3}((-1.0, 1.0, 1.0))
left = Vec{3}((-75.,-25.,-2.))
right = Vec{3}((75.,25.,2.))
model = LandauModel(α, G, (50, 50, 2), left, right, element_potential)

vtk_save("landauorig", model)
@time minimize!(model)
vtk_save("landaufinal", model)

# as we can see this runs very quickly even for relatively large gridsizes.
# The key to get high performance like this is to minimize the allocations inside the threaded loops,
# ideally to 0.


text:
# Implements the Luenberger productivity growth indicator
immutable Luenberger{S<:AbstractDataEnvelopment,T<:RS}
  D0::DDF
  D1::DDF
  K::Int64

  function Luenberger{S,T}(X0::Array,Y0::Array,gx0::Array,gy0::Array,X1::Array,Y1::Array,gx1::Array,gy1::Array) where {S<:AbstractDataEnvelopment,T<:RS}
    if(size(X0) != size(X1) || size(Y0) != size(Y1))
      error("The number of DMUs must be the same in both periods!")
    end
    new(DDF{S,T}(X0,Y0,gx0,gy0),DDF{S,T}(X1,Y1,gx1,gy1),size(X0,1))
  end

  function Luenberger{S,T}(DDF0::DDF,DDF1::DDF) where {S<:AbstractDataEnvelopment,T<:RS}
    new(DDF0,DDF1)
  end
end

function (L::Luenberger)()
  val = Array{Float64}(L.K)
  for k=1:L.K
    val[k] = -geteff(L.D0(L.D1[k]...)) + geteff(L.D1(L.D0[k]...))
  end
  return (val + geteff(L.D0()) - geteff(L.D1()))./2
end

function TEI(L::Luenberger)
  return geteff(L.D0()) - geteff(L.D1())
end

function TC(L::Luenberger)
  TC0 = Array{Float64}(L.K)
  TC1 = Array{Float64}(L.K)
  for k=1:L.K
      # Compute technical change (TC0) from period 0 to period 1 using observations at time 0
      TC0[k] = geteff(L.D1(L.D0[k]...)) - geteff(L.D0(L.D0[k]...))
      # Compute technical change (TC1) from period 0 to period 1 using observations at time 1
      TC1[k] = geteff(L.D1(L.D1[k]...)) - geteff(L.D0(L.D1[k]...))
  end

  # Technical change (TC) is an arithmetic average of TC0 and TC1
  return (TC0+TC1)./2
end

max_stars_repo_name: kepiej/NEPA


text:
module TestModify

include("preamble.jl")
using BangBang.Extras: modify!!

function modify!(f, d, k)
    d′, v = modify!!(f, d, k)
    @test d′ === d
    return v
end

@testset "modify!(f, ::$constructor, key)" for constructor in [
    Dict,
    # SingletonDict,
]
    @testset "update" begin
        dict = constructor(Dict("a" => 1))

        @test modify!(dict, "a") do old
            Some(something(old, 0) + 1)
        end == Some(2)

        @test Dict(dict) == Dict("a" => 2)
    end

    @testset "insert" begin
        dict = constructor(Dict())

        @test modify!(dict, "a") do old
            Some(something(old, 0) + 1)
        end == Some(1)

        @test Dict(dict) == Dict("a" => 1)
    end

    @testset "delete" begin
        dict = constructor(Dict("a" => 1))
        @test modify!(_ -> nothing, dict, "a") === nothing
        @test Dict(dict) == Dict()
    end

    @testset "no-op" begin
        dict = constructor(Dict("a" => 1))
        @test modify!(_ -> nothing, dict, "b") === nothing
        @test Dict(dict) == Dict("a" => 1)

        dict = constructor(Dict("a" => 1))
        @test modify!(_ -> nothing, dict, :b) === nothing
        @test Dict(dict) == Dict("a" => 1)

        dict = constructor(Dict("a" => 1))
        @test modify!(identity, dict, "a") === Some(1)
        @test Dict(dict) == Dict("a" => 1)

        dict = constructor(Dict(:a => 1))
        @test modify!(identity, dict, :a) === Some(1)
        @test Dict(dict) == Dict(:a => 1)
    end

    @testset "mutation inside `f`" begin
        dict = constructor(Dict())

        @test modify!(dict, "a") do old
            dict["a"] = 0
            Some(something(old, 0) + 1)
        end == Some(1)

        @test Dict(dict) == Dict("a" => 1)
    end
end

@testset "modify!!(f, ::Dict, key)" begin
    @testset "update" begin
        dict = Dict("a" => 1)
        backup = deepcopy(dict)
        dict′, ret = modify!!(dict, "a") do old
            Some(something(old, 0) + 0.5)
        end
        @test ret === Some(1.5)
        @test dict′ == Dict("a" => 1.5)
        @test dict == backup
    end

    @testset "insert" begin
        dict = Dict{Union{},Union{}}()
        backup = deepcopy(dict)
        dict′, ret = modify!!(dict, "a") do old
            Some(something(old, 0) + 1)
        end
        @test ret == Some(1)
        @test dict′ == Dict("a" => 1)
        @test dict == backup
    end

    @testset "no-op" begin
        dict = Dict("a" => 1)
        backup = deepcopy(dict)
        dict′, ret = modify!!(_ -> nothing, dict, :b)
        @test ret === nothing
        @test dict′ === dict == backup
    end

    @testset "mutation inside `f` while widening value type" begin
        dict = Dict{String,Int}()
        dict′, ret = modify!!(dict, "a") do old
            dict["a"] = 0
            Some(something(old, 0) + 0.5)
        end
        @test ret === Some(0.5)
        @test dict′ == Dict("a" => 0.5)
        @test dict == Dict("a" => 0)
    end

    @testset "mutation inside `f` while widening key type" begin
        dict = Dict{String,Int}()
        dict′, ret = modify!!(dict, :b) do old
            dict["a"] = 0
            Some(something(old, 0) + 0.5)
        end
        @test ret === Some(0.5)
        @test dict′ == Dict("a" => 0, :b => 0.5)
        @test dict == Dict("a" => 0)
    end
end

@testset "modify!!(f, ::NamedTuple, key)" begin
    @testset "update" begin
        @test modify!!((a = 1,), :a) do old
            Some(something(old, 0) + 1)
        end === ((a = 2,), Some(2))
    end

    @testset "insert" begin
        @test modify!!(NamedTuple(), :a) do old
            Some(something(old, 0) + 1)
        end === ((a = 1,), Some(1))
    end

    @testset "delete" begin
        @test modify!!((a = 1,), :a) do old
            nothing
        end === (NamedTuple(), nothing)
    end

    @testset "no-op" begin
        @test modify!!((a = 1,), :b) do old
            nothing
        end === ((a = 1,), nothing)
    end
end

end  # module


text:
# Lazy structure for an array whose real and fourier part can be accessed by A.real and A.fourier
# For normalization conventions, see PlaneWaveBasis.jl

import Base.getproperty, Base.setproperty!, Base.propertynames
import Base.*, Base.+, Base./
import Base.eltype

## TODO these functions are minimally optimized for the moment
## TODO this always allocates both real and fourier parts, which might be a bit wasteful
## TODO implement broadcasting. This seems to be non-trivial, so I just implement a few simple methods
## TODO implement compressed representation in Fourier space of real arrays

@enum RFA_state RFA_only_real RFA_only_fourier RFA_both

# TODO make this concrete by parametrizing on the type of the arrays
"""
A structure to facilitate manipulations of an array of type T in both real
and fourier space. Create with `from_real` or `from_fourier`, and access
with `A.real` and `A.fourier`.
"""
mutable struct RealFourierArray{Treal <: Real, T <: Union{Treal, Complex{Treal}}}
# Treal is the underlying real type
# T is the type of the array in real space
    basis::PlaneWaveBasis{Treal}
    _real::AbstractArray{T, 3}
    _fourier::AbstractArray{Complex{Treal}, 3}
    _state::RFA_state
end
# Type of the real part
Base.eltype(A::RealFourierArray{Treal, T}) where {Treal, T} = T

# Constructors
RealFourierArray(basis, real, fourier) = RealFourierArray(basis, real, fourier, 2)
function RealFourierArray(basis::PlaneWaveBasis{T}; iscomplex=false) where {T}
    from_real(basis, zeros(iscomplex ? complex(T) : T, basis.fft_size...))
end

function from_real(basis, real_part::AbstractArray{T}) where {T <: Number}
    RealFourierArray{real(T), T}(basis, real_part, similar(real_part, complex(T)), RFA_only_real)
end
function from_fourier(basis, fourier_part::AbstractArray{T}; assume_real=false) where {T <: Complex}
    if assume_real
        RealFourierArray{real(T), real(T)}(basis, similar(fourier_part, real(T)), fourier_part, RFA_only_fourier)
    else
        RealFourierArray{real(T), T}(basis, similar(fourier_part), fourier_part, RFA_only_fourier)
    end
end

check_real(A::RealFourierArray{Treal, T}) where {Treal, T <: Real} = nothing
check_real(A::RealFourierArray) = check_real(A.real)

function Base.propertynames(array::RealFourierArray, private=false)
    ret = [:basis, :real, :fourier]
    private ? append!(ret, fieldnames(RealFourierArray)) : ret
end

function Base.getproperty(A::RealFourierArray, x::Symbol)
    if x == :real
        if A._state == RFA_only_fourier
            r = G_to_r(A.basis, A._fourier)
            if eltype(A) <: Real
                A._real .= real(r)
            else
                A._real .= r
            end
            setfield!(A, :_state, RFA_both)
        end
        return A._real
    elseif x == :fourier
        if A._state == RFA_only_real
            r_to_G!(A._fourier, A.basis, complex.(A._real))
            setfield!(A, :_state, RFA_both)
        end
        return A._fourier
    else
        getfield(A, x)
    end
end
function Base.setproperty!(A::RealFourierArray, x)
    error("RealFourierArray is intended to be read-only." *
          "(This can be bypassed by `setfield!` if you really want to).")
end

# Algebraic operations
function Base.:+(A::RealFourierArray, B::RealFourierArray) where T
    @assert A.basis == B.basis
    from_real(A.basis, A.real + B.real)
end
function Base.:-(A::RealFourierArray, B::RealFourierArray) where T
    @assert A.basis == B.basis
    from_real(A.basis, A.real - B.real)
end
function Base.:*(α::Number, A::RealFourierArray)
    from_real(A.basis, α * A.real)
end
Base.:*(A::RealFourierArray, α::Number) = α * A
Base.:/(A::RealFourierArray, α::Number) = A * inv(α)


text:
using Pkg

Pkg.update()


text:
global 模具cache=Dict{Symbol,Vector{UInt8}}()
function load_模具(s::Symbol;tarfile="data/shape/$s.shape")
	if haskey(模具cache,s)
		return @inbounds 模具cache[s]
	else
		io=open(tarfile,"r")
		v=Vector{UInt8}()
		sizehint!(v,256)
		readbytes!(io,v,256)
		close(io)
		模具cache[s]=v
		return v
	end
end
function fill_模具(v::Vector{UInt8},select::NTuple{4,ARGB32})
	mat=Matrix{ARGB32}(undef,32,32)
	for id in 0:255
		uint=v[id+1]
		i=id>>3+1
		j1=(id&7)<<2+1
		mat[i,j1]=select[uint&3+1]
		mat[i,j1+1]=select[uint>>2&3+1]
		mat[i,j1+2]=select[uint>>4&3+1]
		mat[i,j1+3]=select[uint>>6+1]
	end
	return mat
end


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
# min

Base.min(v::Real, n::ANode) = min(n, v)
Base.min(n::ANode, v::Real) = ANode(:min_c, (n,), size(n), v)

function Base.min(a::ANode, b::ANode)
  @assert size(a) == size(b)
  ANode(:min, (a,b), size(a))
end

function _min{T,N}(n::Csize_t, y::Array{T,N}, a::Array{T,N}, b::T)
  @assert length(a) == length(y) == n
  @simd for i in 1:n
    @inbounds y[i] = min(a[i], b)
  end
end

function _min{T,N}(n::Csize_t, y::Array{T,N}, a::Array{T,N}, b::Array{T,N})
  @assert length(a) == length(b) == length(y) == n
  @simd for i in 1:n
    @inbounds y[i] = min(a[i], b[i])
  end
end

@cuda_gsl _min{T<:Real}(n::Csize_t, y::&T, a::&T, b::T) "y[i] = min(a[i], b)"
@cuda_gsl _min{T<:Real}(n::Csize_t, y::&T, a::&T, b::&T) "y[i] = min(a[i], b[i])"

do_forward!(n::ANode{:min_c}, y, a) = _min(Csize_t(length(y)), y, a, convert(eltype(y), n.arg))
do_forward!(n::ANode{:min}, y, a, b) = _min(Csize_t(length(y)), y, a, b)

gradient_node(n::ANode{:min_c}, wrt::ANode, b::ANode) = b.*(n .!= n.arg)
gradient_node(n::ANode{:min}, wrt::ANode, b::ANode) = b.*(n .== wrt)

# max
Base.max(v::Real, n::ANode) = max(n, v)
Base.max(n::ANode, v::Real) = ANode(:max_c, (n,), size(n), v)

function Base.max(a::ANode, b::ANode)
  @assert size(a) == size(b)
  ANode(:max, (a,b), size(a))
end

function _max{T,N}(n::Csize_t, y::Array{T,N}, a::Array{T,N}, b::T)
  @assert length(a) == length(y) == n
  @simd for i in 1:n
    @inbounds y[i] = max(a[i], b)
  end
end

function _max{T,N}(n::Csize_t, y::Array{T,N}, a::Array{T,N}, b::Array{T,N})
  @assert length(a) == length(b) == length(y) == n
  @simd for i in 1:n
    @inbounds y[i] = max(a[i], b[i])
  end
end

@cuda_gsl _max{T<:Real}(n::Csize_t, y::&T, a::&T, b::T) "y[i] = max(a[i], b)"
@cuda_gsl _max{T<:Real}(n::Csize_t, y::&T, a::&T, b::&T) "y[i] = max(a[i], b[i])"

do_forward!(n::ANode{:max_c}, y, a) = _max(Csize_t(length(y)), y, a, convert(eltype(y), n.arg))
do_forward!(n::ANode{:max}, y, a, b) = _max(Csize_t(length(y)), y, a, b)

gradient_node(n::ANode{:max_c}, wrt::ANode, b::ANode) = b.*(n .!= n.arg)
gradient_node(n::ANode{:max}, wrt::ANode, b::ANode) = b.*(n .== wrt)

# min (const)

mutates(::ANode{:min!}) = true

min!(v::Real, n::ANode) = min!(n, v)
min!(n::ANode, v::Real) = ANode(:min!, (n,), size(n), v)

function _min{T}(n::Csize_t, x::Array{T}, v::T)
  @assert length(x) == n
  for i in 1:n
    @inbounds if x[i] > v x[i] = v end
  end
end

@cuda_gsl _min{T<:Real}(n::Csize_t, x::&T, v::T) "if (x[i] > v) x[i] = v"

do_forward!(n::ANode{:min!}, x) = _min(Csize_t(length(x)), x, convert(eltype(x), arg(n)))
gradient_node(n::ANode{:min!}, wrt::ANode, b::ANode) = b.*(n .!= n.arg)

# max (const)

mutates(::ANode{:max!}) = true

max!(v::Real, n::ANode) = max!(n, v)
max!(n::ANode, v::Real) = ANode(:max!, (n,), size(n), v)

function _max{T}(n::Csize_t, x::Array{T}, v::T)
  @assert length(x) == n
  for i in 1:n
    @inbounds if x[i] < v x[i] = v end
  end
end

@cuda_gsl _max{T<:Real}(n::Csize_t, x::&T, v::T) "if (x[i] < v) x[i] = v"

do_forward!(n::ANode{:max!}, x) = _max(Csize_t(length(x)), x, convert(eltype(x), arg(n)))
gradient_node(n::ANode{:max!}, wrt::ANode, b::ANode) = b.*(n .!= n.arg)

# clip

mutates(::ANode{:clip!}) = true

function clip{T}(n::Csize_t, x::Array{T}, v::T)
  @assert length(x) == n
  for i in 1:n
    @inbounds if abs(x[i]) > v x[i] = copysign(v, x[i]) end
  end
end

@cuda_gsl clip{T<:Real}(n::Csize_t, x::&T, v::T) "if (fabs(x[i]) > v) x[i] = copysignf(v, x[i])"

clip!(n::ANode, v) = ANode(:clip!, (x,), size(x), v)
do_forward!(n::ANode{:clip!}, x) = clip(x, convert(eltype(x), arg(n)))


text:
using Equate
using Statistics

@testset "Mean" begin
    res = Mean(sgtabX, sgtabY)
    popXonY = ExpandTable(res.table.lYx, ACTmath.xcount)
    @test round(mean(popXonY); digits = 4) == 18.9798
end

@testset "Linear" begin
    res = Linear(sgtabX, sgtabY)
    popXonY = ExpandTable(res.table.lYx, ACTmath.xcount)
    @test round(mean(popXonY); digits = 4) == 18.9798
    @test round(std(popXonY; corrected = false); digits = 4) == 8.9394 # 8.9393 in KB p.50
end

@testset "Equipercentile (lower)" begin
    res = Equipercentile(sgtabX, sgtabY)
    popXonY = ExpandTable(res.table.eYx, ACTmath.xcount)
    @test round(mean(popXonY); digits = 4) == 18.9799
    @test round(std(popXonY; corrected = false); digits = 4) == 8.9352
end

@testset "Percentile Rank Function" begin
    P̂ = round2.(Equate.PRF.(0:1:40, Ref(sgtabX)); digits = 2)
    @test sum(KBp48 .== P̂) == 41
end


text:
"""
    RungeKuttaSolution <: AbstractRungeKuttaSolution

A composite type for an [`AbstractRungeKuttaSolution`](@ref).

# Constructors
```julia
RungeKuttaSolution(u, t[, k])
RungeKuttaSolution(problem, solver; savestages = false)
```

# Arguments
- `u          :: AbstractVector`              : numerical solution.
- `t          :: AbstractVector`              : time grid.
- `k          :: AbstractVector`              : vector of stages.
- `problem    :: AbstractInitialValueProblem` : initial value problem.
- `solver     :: AbstractRungeKuttaSolver`    : Runge-Kutta solver.
- `savestages :: Bool`                        : flags when to save all stages into `k`.

# Functions
- [`getindex` ](@ref) : get value(s) and time.
- [`lastindex`](@ref) : last index.
- [`length`   ](@ref) : number of time steps.
- [`setindex!`](@ref) : set value(s) and time.
- [`show`     ](@ref) : shows name and contents.
- [`size`     ](@ref) : number of variables and time steps.
- [`summary`  ](@ref) : shows name.
"""
mutable struct RungeKuttaSolution{u_T, t_T, k_T} <: AbstractRungeKuttaSolution
    u::u_T
    t::t_T
    k::k_T
end

RungeKuttaSolution(u, t) = RungeKuttaSolution(u, t, nothing)

function RungeKuttaSolution(
    problem::AbstractInitialValueProblem,
    solver::AbstractRungeKuttaSolver;
    savestages::Bool = false
)
    @↓ u0, (t0, tN) ← tspan = problem
    @↓ h = solver.stepsize
    N = round(Int, (tN - t0) / h) + 1
    # Initialise u
    u0_T = eltype(u0)
    L = length(u0)
    u = Vector{u0_T}(undef, N, L)
    u[1] = u0
    # Initialise t
    t0_T = typeof(t0)
    t = Vector{t0_T}(undef, N)
    t[1] = t0
    # Initialise k
    if savestages
        @↓ s = solver.tableau
        k = Vector{u0_T}(undef, N, s, L)
        return RungeKuttaSolution(u, t, k)
    else
        return RungeKuttaSolution(u, t)
    end
end

############################################################################################
#                                         FUNCTIONS                                        #
############################################################################################

"""
    length(solution::AbstractRungeKuttaSolution)

returns the number of time steps of `solution`.
"""
Base.length(solution::AbstractRungeKuttaSolution) = length(solution.t)

"""
    size(solution::AbstractRungeKuttaSolution)

returns a tuple containing the number of variables and of time steps of `solution`.
"""
Base.size(solution::AbstractRungeKuttaSolution) = (length(solution.u[1]), length(solution.t))

"""
    getindex(solution::AbstractRungeKuttaSolution, i::Integer)

returns a [`RungeKuttaSolution`](@ref) containing the fields of `solution` indexed at `i`.
"""
function Base.getindex(solution::AbstractRungeKuttaSolution, i::Integer)
    @↓ u, t, k = solution
    if k isa Nothing
        return RungeKuttaSolution(u[i], t[i])
    else
        return RungeKuttaSolution(u[i], t[i], k[i])
    end
end

"""
    setindex!(solution::AbstractRungeKuttaSolution, tuple::Tuple, i::Integer)

stores the values of `tuple` into the fields of `solution` indexed at `i`.
"""
function Base.setindex!(solution::AbstractRungeKuttaSolution, tuple::Tuple, i::Integer)
    @↓ u, t, k = solution
    u[i] = tuple[1]
    t[i] = tuple[2]
    if length(tuple) > 2 && !(k isa Nothing)
        k[i] = tuple[3]
    end
    return solution
end

"""
    lastindex(solution::AbstractRungeKuttaSolution)

returns the last index of `solution`.
"""
Base.lastindex(solution::AbstractRungeKuttaSolution) = lastindex(solution.t)

############################################################################################
#                                         PRINTING                                         #
############################################################################################

"""
    show(io::IO, solution::AbstractRungeKuttaSolution)

prints a full description of `solution` and its contents to a stream `io`.
"""
Base.show(io::IO, solution::AbstractRungeKuttaSolution) = NSDEBase._show(io, solution)

"""
    summary(io::IO, solution::AbstractRungeKuttaSolution)

prints a brief description of `solution` to a stream `io`.
"""
Base.summary(io::IO, solution::AbstractRungeKuttaSolution) = NSDEBase._summary(io, solution)

############################################################################################
#                                          METHODS                                         #
############################################################################################

function extract(solution::AbstractRungeKuttaSolution, i)
    (L, N) = size(solution)
    @↓ u, t = solution
    return RungeKuttaSolution([u[n][i] for n = 1:N], t)
end

# Make solution callable
function findclosest(target, vector)
    N = length(vector)
    if target ≤ vector[1]
        return (1, vector[1])
    elseif target ≥ vector[end]
        return (N, vector[end])
    end
    i = 0; j = N; n = 0
    while i < j
        n = (i + j) ÷ 2
        if target == vector[n]
            return vector[n]
        elseif target < vector[n]
            if (n > 1) && (target > vector[n - 1])
                if target - vector[n - 1] ≥ vector[n] - target
                    return (n, vector[n])
                else
                    return (n - 1, vector[n - 1])
                end
            end
            j = n
        else
            if (n < N) && (target < vector[n + 1])
                if target - vector[n] ≥ vector[n + 1] - target
                    return (n + 1, vector[n + 1])
                else
                    return (n, vector[n])
                end
            end
            i = n + 1
        end
    end
    return n, vector[n]
end

function (solution::AbstractRungeKuttaSolution)(tₚ::Real)
    @↓ u, t = solution
    (n, tₙ) = findclosest(tₚ, t)
    return tₙ, u[n]
end

function (solution::AbstractRungeKuttaSolution)(tspan::Tuple{Real, Real})
    @↓ u, t = solution
    n₁, _ = findclosest(tspan[1], t)
    n₂, _ = findclosest(tspan[2], t)
    return RungeKuttaSolution(u[n₁:n₂], t[n₁:n₂])
end

(solution::AbstractRungeKuttaSolution)(t₁::Real, t₂::Real) = solution((t₁, t₂))


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

struct ConvexPolygon{T,N,M} <: ParametricSurface{T,N}
    plane::Plane{T,N}
    origin::SVector{N,T}
    vertices::SMatrix{N,M,T}
end

# uv(surface::ParametricSurface{T,N}, p::SVector{N,T}) -> SVector{2,T}
# uvrange(surface::ParametricSurface{T,N}) -> Tuple{Tuple{T,T},Tuple{T,T}}
# point(surface::ParametricSurface{T,N}, u::T, v::T) -> SVector{N,T}
# partials(surface::ParametricSurface{T,N}, u::T, v::T) -> Tuple{SVector{N,T}, SVector{N,T}}
# normal(surface::ParametricSurface{T,N}, u::T, v::T) -> SVector{N,T}
# inside(surface::ParametricSurface{T,N}, p: :SVector{N,T}) -> Bool
# onsurface(surface::ParametricSurface{T,N}, p::SVector{N,T}) -> Bool
# surfaceintersection(surface::ParametricSurface{T,N}, AbstractRay::Ray{T,N}) -> Union{EmptyInterval{T},Interval{T},DisjointUnion{T}}
# interface(surface::ParametricSurface{T,N}) -> OpticalInterface{T}