Skip to content

Instantly share code, notes, and snippets.

Avatar
🥝
Focusing

Yingbo Ma YingboMa

🥝
Focusing
View GitHub Profile
View periodic_qr.jl
using LinearAlgebra, BandedMatrices
n = 3000
A = Matrix(SymTridiagonal(Symmetric(rand(n, n)))); A[end, 1] = A[1, end] = 1;
b = rand(size(A, 1))
function periodic_givens_rotation!(A::AbstractMatrix{T}) where T
m, n = size(A)
gs = Vector{LinearAlgebra.Givens{T}}(undef, 2n-2)
jj = 0
View clenshaw_vs_gauss.jl
using Plots, FastGaussQuadrature, FastTransforms
fejer1(N) = FastTransforms.fejernodes1(Float64, N), FastTransforms.fejerweights1(FastTransforms.FastTransforms.chebyshevmoments1(Float64, N))
fejer2(N) = FastTransforms.fejernodes2(Float64, N), FastTransforms.fejerweights2(FastTransforms.FastTransforms.chebyshevmoments2(Float64, N))
clenshawcurtis(N) = FastTransforms.clenshawcurtisnodes(Float64, N), FastTransforms.clenshawcurtisweights(FastTransforms.FastTransforms.chebyshevmoments1(Float64, N))
x = range(-1.2, stop = 1.2, length = 1000)
y = range(-0.4, stop = 0.4, length = 1000÷3)
r(z, x, w) = sum(k->w[k]/(z - x[k]), eachindex(x))
ϕ(z) = log((z+1)/(z-1))
quadrature_contour(xs, w; kwargs...) = contour(x, y, (x, y) -> abs(r(x + y*im, xs, w) - ϕ(x + y*im)); levels=map(i->1/10^i, 0:10), aspect_ratio=1, colorbar=false, xlimits=(-1.5, 1.5), ylimits=(-0.5, 0.5), kwargs...)
plts = [quadrature_contour(clenshawcurtis(64)...; title="Clenshaw-Curtis"),
View stability_contour.jl
using DiffEqDevTools, BifurcationKit, Setfield, Plots
mutable struct TermCache
start
sign::Bool
count::Int
end
function stability_curve(tab)
function curve((x, ), (y, ))
View teach_forwarddiff.jl
using ForwardDiff
goo((x, y, z),) = [x^2*z, x*y*z, abs(z)-y]
foo((x, y, z),) = [x^2*z, x*y*z, abs(z)-y]
function foo(u::Vector{ForwardDiff.Dual{T,V,P}}) where {T,V,P}
# unpack: AoS -> SoA
vs = ForwardDiff.value.(u)
# you can play with the dimension here, sometimes it makes sense to transpose
ps = mapreduce(ForwardDiff.partials, hcat, u)
# get f(vs)
val = foo(vs)
View linear_solve_design.jl
#=
Replace calcJ, calcW, linsolve, do_newJW to [u, p, t, newJ, newW, tol, error_estimate]
foo!(z, integrator, alg, nlsolver) -> z .= W\z
foo(z, integrator, alg, nlsolver) -> W\z
foo!(z, integrators, alg, nlsolver)
- User control inverse of W(_t)
- user defined factorization (object with ldiv! defined) (Info: u, p, t, newW,
tol, error_estimate)
@YingboMa
YingboMa / radau.jl
Created Apr 3, 2020
Radau IIA transformed coefficients
View radau.jl
A = [11//45-7sqrt(6)/360 37//225-169sqrt(6)/1800 -2//225+sqrt(6)/75
37//225+169sqrt(6)/1800 11//45+7sqrt(6)/360 -2//225-sqrt(6)/75
4//9-sqrt(6)/36 4//9+sqrt(6)/36 1//9]
T = Float64
T11 = convert(T, 9.1232394870892942792e-02)
T12 = convert(T, -0.14125529502095420843e0)
T13 = convert(T, -3.0029194105147424492e-02)
T21 = convert(T, 0.24171793270710701896e0)
T22 = convert(T, 0.20412935229379993199e0)
View fd_typed.jl
│ ─ %-1 = invoke perform_step!(::OrdinaryDiffEq.ODEIntegrator{Tsit5,true,Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float64,4},1},Float64,DiffEqBase.NullParameters,Float
64,Float64,Float64,Array{Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float64,4},1},1},ODESolution{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float64,4},2,Array{
Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float64,4},1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float
64,4},1},1},1},ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{var"#25#26",Float64},Float64,4},1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(lorenz),L
inearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple
{}}},DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFuncti
View fd_llvm.ll
define void @"julia_perform_step!_19254"(%jl_value_t addrspace(10)* nonnull align 8 dereferenceable(304), %jl_value_t addrspace(10)* nonnull align 8 dereferenceable(544), i8) {
top:
%3 = addrspacecast %jl_value_t addrspace(10)* %0 to %jl_value_t addrspace(11)*
%4 = bitcast %jl_value_t addrspace(11)* %3 to i8 addrspace(11)*
%5 = getelementptr inbounds i8, i8 addrspace(11)* %4, i64 32
%6 = bitcast i8 addrspace(11)* %5 to double addrspace(11)*
%7 = load double, double addrspace(11)* %6, align 8
%8 = getelementptr inbounds i8, i8 addrspace(11)* %4, i64 48
%9 = bitcast i8 addrspace(11)* %8 to %jl_value_t addrspace(10)* addrspace(11)*
%10 = load %jl_value_t addrspace(10)*, %jl_value_t addrspace(10)* addrspace(11)* %9, align 8
@YingboMa
YingboMa / specialize_vararg.jl
Created Jan 12, 2020
Got from Mason Protter
View specialize_vararg.jl
# by Mason Protter
using MacroTools: MacroTools, splitdef, combinedef, @capture
macro specialize_vararg(n::Int, fdef)
d = splitdef(fdef)
args = d[:args][end]
@assert d[:args][end] isa Expr && d[:args][end].head == Symbol("...")
args_symbol = d[:args][end].args[]
fdefs = Expr(:block)
for i in 1:n-1
di = deepcopy(d)
View D_sicm.jl
using ModelingToolkit
using ModelingToolkit: Differential, expand_derivatives, Expression, Operation, simplify_constants
function D(f)
(args...) -> begin
syms = map(_->Variable(gensym())(), args)
ex = f(syms...)
ds = map(s->Differential(s), syms)
ops = [expand_derivatives(ds[i](ex)) for i in eachindex(syms)]
rules = Dict( collect(zip(syms, args)) )