diff --git a/src/ACME.jl b/src/ACME.jl index 85de2395..ffbec765 100644 --- a/src/ACME.jl +++ b/src/ACME.jl @@ -110,7 +110,7 @@ include("elements.jl") include("circuit.jl") -mutable struct DiscreteModel{Solvers} +mutable struct DiscreteModel{Solvers<:Tuple{Vararg{NonlinearSolver}}} a::Matrix{Float64} b::Matrix{Float64} c::Matrix{Float64} @@ -131,18 +131,14 @@ mutable struct DiscreteModel{Solvers} solvers::Solvers x::Vector{Float64} - function DiscreteModel(mats::Dict{Symbol}, nonlinear_eq_funcs::Vector, + function DiscreteModel(mats, nonlinear_eq_funcs::Vector, solvers::Solvers) where {Solvers} - model = new{Solvers}() - - for mat in (:a, :b, :c, :pexps, :dqs, :eqs, :fqprevs, :fqs, :dy, :ey, :fy, :x0, :q0s, :y0) - setfield!(model, mat, convert(fieldtype(typeof(model), mat), mats[mat])) - end - - model.nonlinear_eq_funcs = nonlinear_eq_funcs - model.solvers = solvers - model.x = zeros(nx(model)) - return model + return new{Solvers}( + mats[:a], mats[:b], mats[:c], mats[:x0], + mats[:pexps], mats[:dqs], mats[:eqs], mats[:fqprevs], mats[:fqs], mats[:q0s], + mats[:dy], mats[:ey], mats[:fy], mats[:y0], + nonlinear_eq_funcs, solvers, zeros(length(mats[:x0])) + ) end end @@ -159,12 +155,11 @@ function DiscreteModel(circ::Circuit, t::Real, ::Type{Solver}=HomotopySolver{Cac end model_nns = Int[sum(nns[nles]) for nles in nl_elems] - model_qidxs = [vcat(consecranges(nqs)[nles]...) for nles in nl_elems] - split_nl_model_matrices!(mats, model_qidxs, model_nns) + model_qidxs = [reduce(vcat, consecranges(nqs)[nles]) for nles in nl_elems] + mats = merge(mats, split_nl_model_matrices(mats, model_qidxs, model_nns)) - reduce_pdims!(mats) + mats = reduce_pdims!(mats) - model_nps = size.(mats[:dqs], 1) model_nqs = size.(mats[:pexps], 1) @assert nn(circ) == sum(model_nns) @@ -173,7 +168,7 @@ function DiscreteModel(circ::Circuit, t::Real, ::Type{Solver}=HomotopySolver{Cac fqs = Matrix{Float64}.(mats[:fqs]) fqprev_fulls = Matrix{Float64}.(mats[:fqprev_fulls]) - model_nonlinear_eq_funcs = [ + model_nonlinear_eq_funcs = Function[ let q = zeros(nq), circ_nl_func = nonlinear_eq_func(circ, nles) @inline function(res, J, pfull, Jq, fq, z) #copyto!(q, pfull + fq * z) @@ -188,27 +183,30 @@ function DiscreteModel(circ::Circuit, t::Real, ::Type{Solver}=HomotopySolver{Cac end end for (nles, nq) in zip(nl_elems, model_nqs)] - nonlinear_eq_funcs = [ + nonlinear_eq_funcs = Function[ @inline function (res, J, scratch, z) nleq(res, J, scratch[1], scratch[2], fq, z) end for (nleq, fq) in zip(model_nonlinear_eq_funcs, fqs)] init_zs = [zeros(nn) for nn in model_nns] for idx in eachindex(nonlinear_eq_funcs) - q = q0s[idx] + fqprev_fulls[idx] * vcat(init_zs...) + q = q0s[idx] + fqprev_fulls[idx] * reduce(vcat, init_zs) init_zs[idx] = initial_solution(nonlinear_eq_funcs[idx], q, model_nns[idx]) end - while any(np -> np == 0, model_nps) - const_idxs = findall(iszero, model_nps) - const_zidxs = vcat(consecranges(model_nns)[const_idxs]...) + while true + const_idxs = findall(iszero, size.(mats[:dqs], 1)) + if isempty(const_idxs) + break + end + const_zidxs = reduce(vcat, consecranges(model_nns)[const_idxs]) varying_zidxs = filter(idx -> !(idx in const_zidxs), 1:sum(model_nns)) for idx in eachindex(mats[:q0s]) - mats[:q0s][idx] += mats[:fqprev_fulls][idx][:,const_zidxs] * vcat(init_zs[const_idxs]...) + mats[:q0s][idx] += mats[:fqprev_fulls][idx][:,const_zidxs] * reduce(vcat, init_zs[const_idxs]) mats[:fqprev_fulls][idx] = mats[:fqprev_fulls][idx][:,varying_zidxs] end - mats[:x0] += mats[:c][:,const_zidxs] * vcat(init_zs[const_idxs]...) - mats[:y0] += mats[:fy][:,const_zidxs] * vcat(init_zs[const_idxs]...) + copyto!(mats[:x0], mats[:x0] + mats[:c][:,const_zidxs] * reduce(vcat, init_zs[const_idxs])) + copyto!(mats[:y0], mats[:y0] + mats[:fy][:,const_zidxs] * reduce(vcat, init_zs[const_idxs])) deleteat!(mats[:q0s], const_idxs) deleteat!(mats[:dq_fulls], const_idxs) deleteat!(mats[:eq_fulls], const_idxs) @@ -220,11 +218,10 @@ function DiscreteModel(circ::Circuit, t::Real, ::Type{Solver}=HomotopySolver{Cac deleteat!(model_nonlinear_eq_funcs, const_idxs) deleteat!(nonlinear_eq_funcs, const_idxs) deleteat!(nl_elems, const_idxs) - mats[:fy] = mats[:fy][:,varying_zidxs] - mats[:c] = mats[:c][:,varying_zidxs] - reduce_pdims!(mats) - model_nps = size.(mats[:dqs], 1) + mats = merge(mats, (fy=mats[:fy][:,varying_zidxs], c=mats[:c][:,varying_zidxs])) + mats = reduce_pdims!(mats) end + model_nps = size.(mats[:dqs], 1) q0s = Array{Float64}.(mats[:q0s]) fqs = Array{Float64}.(mats[:fqs]) @@ -260,50 +257,56 @@ function DiscreteModel(circ::Circuit, t::Real, ::Type{Solver}=HomotopySolver{Cac end function model_matrices(circ::Circuit, t::Rational{BigInt}) - lhs = convert(SparseMatrixCSC{Rational{BigInt},Int}, - [mv(circ) mi(circ) mxd(circ)//t+mx(circ)//2 mq(circ); - blockdiag(topomat(circ)...) spzeros(nb(circ), nx(circ) + nq(circ))]) - rhs = convert(SparseMatrixCSC{Rational{BigInt},Int}, - [u0(circ) mu(circ) mxd(circ)//t-mx(circ)//2; - spzeros(nb(circ), 1+nu(circ)+nx(circ))]) + lhs = Rational{BigInt}[ + mv(circ) mi(circ) mxd(circ)//t+mx(circ)//2 mq(circ); + blockdiag(topomat(circ)...) spzeros(nb(circ), nx(circ) + nq(circ)) + ] + rhs = Rational{BigInt}[ + u0(circ) mu(circ) mxd(circ)//t-mx(circ)//2; + spzeros(Rational{BigInt}, nb(circ), 1+nu(circ)+nx(circ)) + ] x, f = Matrix.(gensolve(lhs, rhs)) - rowsizes = [nb(circ); nb(circ); nx(circ); nq(circ)] - res = Dict{Symbol,Array}(zip([:fv; :fi; :c; :fq], matsplit(f, rowsizes))) + rowsizes = (nb(circ), nb(circ), nx(circ), nq(circ)) + rowranges = consecranges(rowsizes) + fq = f[rowranges[4], :] - nullspace = gensolve(sparse(res[:fq]::Matrix{Rational{BigInt}}), - spzeros(Rational{BigInt}, size(res[:fq],1), 0))[2] + nullspace = gensolve(sparse(fq), spzeros(Rational{BigInt}, size(fq, 1), 0))[2] indeterminates = f * nullspace - if sum(abs2, res[:c] * nullspace) > 1e-20 + if sum(abs2, indeterminates[rowranges[3],:]) > 1e-20 @warn "State update depends on indeterminate quantity" end + while size(nullspace, 2) > 0 - i, j = argmax(abs.(nullspace)).I + i, j = (argmax(abs.(nullspace))::CartesianIndex{2}).I # argmax cannot be inferred prior to Julia 1.7 nullspace = nullspace[[1:i-1; i+1:end], [1:j-1; j+1:end]] f = f[:, [1:i-1; i+1:end]] - for k in [:fv; :fi; :c; :fq] - res[k] = res[k][:, [1:i-1; i+1:end]] - end end - merge!(res, Dict(zip([:v0 :ev :dv; :i0 :ei :di; :x0 :b :a; :q0 :eq_full :dq_full], - matsplit(x, rowsizes, [1; nu(circ); nx(circ)])))) - for v in (:v0, :i0, :x0, :q0) - res[v] = dropdims(res[v], dims=2) - end + f_split = NamedTuple{(:fv, :fi, :c, :fq)}(matsplit(f, rowsizes)) + + x_split = NamedTuple{(:v0, :i0, :x0, :q0, :ev, :ei, :b, :eq_full, :dv, :di, :a, :dq_full)}( + matsplit(x, rowsizes, (1, nu(circ), nx(circ))) + ) + vs = (:v0, :i0, :x0, :q0) + x_split_replace0 = NamedTuple{vs}( + map(let x_split=x_split; v -> dropdims(x_split[v], dims=2); end, vs) + ) p = [pv(circ) pi(circ) px(circ)//2+pxd(circ)//t pq(circ)] if sum(abs2, p * indeterminates) > 1e-20 @warn "Model output depends on indeterminate quantity" end - res[:dy] = p * x[:,2+nu(circ):end] + px(circ)//2-pxd(circ)//t - # p * [dv; di; a; dq_full] + px(circ)//2-pxd(circ)//t - res[:ey] = p * x[:,2:1+nu(circ)] # p * [ev; ei; b; eq_full] - res[:fy] = p * f # p * [fv; fi; c; fq] - res[:y0] = p * vec(x[:,1]) # p * [v0; i0; x0; q0] + y_split = ( + dy = p * x[:,2+nu(circ):end] + px(circ)//2-pxd(circ)//t, + # p * [dv; di; a; dq_full] + px(circ)//2-pxd(circ)//t + ey = p * x[:,2:1+nu(circ)], # p * [ev; ei; b; eq_full] + fy = p * f, # p * [fv; fi; c; fq] + y0 = p * x[:,1], # p * [v0; i0; x0; q0] + ) - return res + return merge(f_split, x_split, x_split_replace0, y_split) end model_matrices(circ::Circuit, t) = model_matrices(circ, Rational{BigInt}(t)) @@ -349,7 +352,7 @@ function nldecompose!(mats, nns, nqs) while !isempty(rem_nles) for sz in 1:length(rem_nles), sub in subsets(collect(rem_nles), sz) nn_sub = sum(nns[sub]) - a_update = tryextract(fq[[sub_ranges[sub]...;],rem_cols], nn_sub) + a_update = tryextract(fq[reduce(vcat, sub_ranges[sub]), rem_cols], nn_sub) if a_update !== nothing fq[:,rem_cols] = fq[:,rem_cols] * a_update a[:,rem_cols] = a[:,rem_cols] * a_update @@ -363,38 +366,49 @@ function nldecompose!(mats, nns, nqs) end end - mats[:c] = mats[:c] * a + copyto!(mats[:c], mats[:c] * a) # mats[:fq] is updated as part of the loop - mats[:fy] = mats[:fy] * a + copyto!(mats[:fy], mats[:fy] * a) return extracted_subs end -function split_nl_model_matrices!(mats, model_qidxs, model_nns) - mats[:dq_fulls] = Matrix[mats[:dq_full][qidxs,:] for qidxs in model_qidxs] - mats[:eq_fulls] = Matrix[mats[:eq_full][qidxs,:] for qidxs in model_qidxs] - let fqsplit = vcat((matsplit(mats[:fq][qidxs,:], [length(qidxs)], model_nns) for qidxs in model_qidxs)...) - mats[:fqs] = Matrix[fqsplit[i,i] for i in 1:length(model_qidxs)] - mats[:fqprev_fulls] = Matrix[[fqsplit[i, 1:i-1]... zeros(eltype(mats[:fq]), length(model_qidxs[i]), sum(model_nns[i:end]))] - for i in 1:length(model_qidxs)] - end - mats[:q0s] = Vector[mats[:q0][qidxs] for qidxs in model_qidxs] +function split_nl_model_matrices(mats, model_qidxs, model_nns) + fqsplit = reduce( + vcat, + matsplit(mats[:fq][qidxs,:], [length(qidxs)], model_nns) for qidxs in model_qidxs; + init=Matrix{typeof(mats[:fq])}(undef, 0, length(model_nns)) + ) + return ( + dq_fulls = [mats[:dq_full][qidxs,:] for qidxs in model_qidxs], + eq_fulls = [mats[:eq_full][qidxs,:] for qidxs in model_qidxs], + fqs = [fqsplit[i,i] for i in 1:length(model_qidxs)], + fqprev_fulls = [ + foldr( + hcat, + fqsplit[i, 1:i-1]; + init = zeros(eltype(mats[:fq]), length(model_qidxs[i]), sum(model_nns[i:end])) + ) + for i in 1:length(model_qidxs) + ], + q0s = [mats[:q0][qidxs] for qidxs in model_qidxs] + ) end -function reduce_pdims!(mats::Dict) +function reduce_pdims!(mats) subcount = length(mats[:dq_fulls]) - mats[:dqs] = Vector{Matrix}(undef, subcount) - mats[:eqs] = Vector{Matrix}(undef, subcount) - mats[:fqprevs] = Vector{Matrix}(undef, subcount) - mats[:pexps] = Vector{Matrix}(undef, subcount) + dqs = Vector{eltype(mats[:dq_fulls])}(undef, subcount) + eqs = Vector{eltype(mats[:eq_fulls])}(undef, subcount) + fqprevs = Vector{eltype(mats[:fqprev_fulls])}(undef, subcount) + pexps = Vector{Matrix{Rational{BigInt}}}(undef, subcount) offset = 0 for idx in 1:subcount # decompose [dq_full eq_full] into pexp*[dq eq] with [dq eq] having minimum # number of rows pexp, dqeq = rank_factorize(sparse([mats[:dq_fulls][idx] mats[:eq_fulls][idx] mats[:fqprev_fulls][idx]])) - mats[:pexps][idx] = pexp - colsizes = [size(mats[m][idx], 2) for m in [:dq_fulls, :eq_fulls, :fqprev_fulls]] - mats[:dqs][idx], mats[:eqs][idx], mats[:fqprevs][idx] = matsplit(dqeq, [size(dqeq, 1)], colsizes) + pexps[idx] = pexp + colsizes = map(m -> size(mats[m][idx], 2)::Int, (:dq_fulls, :eq_fulls, :fqprev_fulls)) + dqs[idx], eqs[idx], fqprevs[idx] = matsplit(dqeq, (size(dqeq, 1),), colsizes) # project pexp onto the orthogonal complement of the column space of Fq fq = mats[:fqs][idx] @@ -403,27 +417,32 @@ function reduce_pdims!(mats::Dict) pexp = pexp - fq*fq_pinv*pexp # if the new pexp has lower rank, update pexp, f = rank_factorize(sparse(pexp)) - if size(pexp, 2) < size(mats[:pexps][idx], 2) + if size(pexp, 2) < size(pexps[idx], 2) cols = offset .+ (1:nn) - mats[:a] = mats[:a] - mats[:c][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:dqs][idx] - mats[:b] = mats[:b] - mats[:c][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:eqs][idx] - mats[:dy] = mats[:dy] - mats[:fy][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:dqs][idx] - mats[:ey] = mats[:ey] - mats[:fy][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:eqs][idx] + copyto!(mats[:a], mats[:a] - mats[:c][:,cols]*fq_pinv*pexps[idx]*dqs[idx]) + copyto!(mats[:b], mats[:b] - mats[:c][:,cols]*fq_pinv*pexps[idx]*eqs[idx]) + copyto!(mats[:dy], mats[:dy] - mats[:fy][:,cols]*fq_pinv*pexps[idx]*dqs[idx]) + copyto!(mats[:ey], mats[:ey] - mats[:fy][:,cols]*fq_pinv*pexps[idx]*eqs[idx]) for idx2 in (idx+1):subcount - mats[:dq_fulls][idx2] = mats[:dq_fulls][idx2] - mats[:fqprev_fulls][idx2][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:dqs][idx] - mats[:eq_fulls][idx2] = mats[:eq_fulls][idx2] - mats[:fqprev_fulls][idx2][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:eqs][idx] - mats[:fqprev_fulls][idx2][:,1:offset] = mats[:fqprev_fulls][idx2][:,1:offset] - mats[:fqprev_fulls][idx2][:,cols]*fq_pinv*mats[:pexps][idx]*mats[:fqprevs][idx][:,1:offset] + q = mats[:fqprev_fulls][idx2][:,cols]*fq_pinv*pexps[idx] + copyto!(mats[:dq_fulls][idx2], mats[:dq_fulls][idx2] - q*dqs[idx]) + copyto!(mats[:eq_fulls][idx2], mats[:eq_fulls][idx2] - q*eqs[idx]) + copyto!( + mats[:fqprev_fulls][idx2][:,1:offset], + mats[:fqprev_fulls][idx2][:,1:offset] - q*fqprevs[idx][:,1:offset] + ) end - mats[:pexps][idx] = pexp - mats[:dqs][idx] = f * mats[:dqs][idx] - mats[:eqs][idx] = f * mats[:eqs][idx] - mats[:fqprevs][idx] = f * mats[:fqprevs][idx] - mats[:dq_fulls][idx] = pexp * mats[:dqs][idx] - mats[:eq_fulls][idx] = pexp * mats[:eqs][idx] - mats[:fqprev_fulls][idx] = pexp * mats[:fqprevs][idx] + pexps[idx] = pexp + dqs[idx] = f * dqs[idx] + eqs[idx] = f * eqs[idx] + fqprevs[idx] = f * fqprevs[idx] + copyto!(mats[:dq_fulls][idx], pexp * dqs[idx]) + copyto!(mats[:eq_fulls][idx], pexp * eqs[idx]) + copyto!(mats[:fqprev_fulls][idx], pexp * fqprevs[idx]) end offset += nn end + return merge(mats, (dqs = dqs, eqs=eqs, fqprevs=fqprevs, pexps=pexps)) end function initial_solution(init_nl_eq_func::Function, q0, nn) @@ -695,7 +714,7 @@ function gensolve(a, b, x, h, thresh=0.1) if m == 0 return x, h end - t = sortperm(vec(mapslices(ait -> count(!iszero, ait), a, dims=2))) # row indexes in ascending order of nnz + t = sortperm([count(!iszero, ait) for ait in eachrow(a)]) # row indexes in ascending order of nnz tol = 3 * max(eps(float(eltype(a))), eps(float(eltype(h)))) * size(a, 2) for i in 1:m ait = a[t[i],:]' # ait is a row of the a matrix @@ -707,7 +726,7 @@ function gensolve(a, b, x, h, thresh=0.1) continue end jat = jnz[nz_abs_vals .≥ thresh*max_abs_val] # cols above threshold - j = jat[argmin(vec(mapslices(hj -> count(!iszero, hj), h[:,jat], dims=1)))] + j = jat[argmin([count(!iszero, hj) for hj in eachcol(h[:,jat])])] q = h[:,j] x = x + convert(typeof(x), q * ((b[t[i],:]' - ait*x) * (1 / (ait*q)))) if size(h)[2] > 1 @@ -727,7 +746,7 @@ function rank_factorize(a::SparseMatrixCSC) nullspace = gensolve(a', spzeros(eltype(a), size(a, 2), 0))[2] c = Matrix{eltype(a)}(I, size(a, 1), size(a, 1)) while size(nullspace, 2) > 0 - i, j = argmax(abs.(nullspace)).I + i, j = (argmax(abs.(nullspace))::CartesianIndex{2}).I # argmax cannot be inferred prior to Julia 1.7 c -= c[:, i] * nullspace[:, j]' / nullspace[i, j] c = c[:, [1:i-1; i+1:end]] nullspace -= nullspace[:, j] * vec(nullspace[i, :])' / nullspace[i, j] @@ -737,9 +756,18 @@ function rank_factorize(a::SparseMatrixCSC) return c, f end -consecranges(lengths) = map((l, e) -> (e-l+1):e, lengths, cumsum(lengths)) +# cumsum(::Tuple) requires Julia 1.5, so roll our own version if needed, named +# _cumsum so not to commit type piracy +_cumsum(x) = cumsum(x) +if !hasmethod(cumsum, Tuple{Tuple}) + _cumsum(x::Tuple) = Base.tail(foldl((r, xi) -> (r..., r[end] + xi), x; init=(false,))) +end + +consecranges(lengths) = map((l, e) -> (e-l+1):e, lengths, _cumsum(lengths)) matsplit(v::AbstractVector, rowsizes) = [v[rs] for rs in consecranges(rowsizes)] +matsplit(m::AbstractMatrix, rowsizes::Tuple, colsizes::Tuple=(size(m,2),)) = + ((map(cs -> map(rs -> m[rs, cs], consecranges(rowsizes)), consecranges(colsizes))...)...,) matsplit(m::AbstractMatrix, rowsizes, colsizes=[size(m,2)]) = [m[rs, cs] for rs in consecranges(rowsizes), cs in consecranges(colsizes)] @@ -768,4 +796,18 @@ precompile(mosfet, (Symbol,)) precompile(opamp, ()) precompile(opamp, (Type{Val{:macak}}, Float64, Float64, Float64)) +precompile(Circuit, ()) +precompile(add!, (Circuit, Element)) +precompile(add!, (Circuit, Symbol, Element)) +for T1 in (Symbol, Tuple{Symbol,Int}, Tuple{Symbol,String}, Tuple{Symbol,Symbol}), + T2 in ( + Symbol, Tuple{Symbol,Int}, Tuple{Symbol,String}, Tuple{Symbol,Symbol}, + Vararg{Symbol}, Vararg{Tuple{Symbol,Int}}, Vararg{Tuple{Symbol,String}}, Vararg{Tuple{Symbol,Symbol}}, + Vararg{Union{Symbol,Tuple{Symbol,Any}}}, + ) + precompile(connect!, (Circuit, T1, T2)) +end + +precompile(DiscreteModel, (Circuit, Rational{Int})) + end # module diff --git a/src/circuit.jl b/src/circuit.jl index f13a67c0..3b36dcf4 100644 --- a/src/circuit.jl +++ b/src/circuit.jl @@ -25,7 +25,7 @@ struct Circuit elements::OrderedDict{Symbol, Element} nets :: Vector{Net} net_names :: Dict{Symbol, Net} - Circuit() = new(OrderedDict{Symbol, Element}(), [], Dict{Symbol, Net}()) + Circuit() = new(OrderedDict{Symbol, Element}(), Net[], Dict{Symbol, Net}()) end elements(c::Circuit) = values(c.elements) @@ -46,7 +46,7 @@ for mat in [:mv; :mi; :mx; :mxd; :mq; :mu; :pv; :pi; :px; :pxd; :pq] )::SparseMatrixCSC{Rational{BigInt},Int} end -u0(c::Circuit) = vcat((elem.u0 for elem in elements(c))...) +u0(c::Circuit) = reduce(vcat, elem.u0 for elem in elements(c); init=spzeros(Real, 0, 1)) function incidence(c::Circuit) i = sizehint!(Int[], 2nb(c)) @@ -153,7 +153,7 @@ end netfor!(c::Circuit, p::Tuple{Symbol,Any}) = netfor!(c, (p[1], Symbol(p[2]))) function netfor!(c::Circuit, name::Symbol) - haskey(c.net_names, name) || push!(c.nets, get!(c.net_names, name, [])) + haskey(c.net_names, name) || push!(c.nets, get!(c.net_names, name, Net[])) c.net_names[name] end @@ -179,10 +179,12 @@ connect!(circ, (:src, -), (:r, 2), :gnd) # connect to gnd net ``` """ function connect!(c::Circuit, pins::Union{Symbol,Tuple{Symbol,Any}}...) - nets = unique([netfor!(c, pin) for pin in pins]) + nets = unique([netfor!(c, pin)::Net for pin in pins]) for net in nets[2:end] append!(nets[1], net) - deleteat!(c.nets, findfirst(n -> n == net, c.nets)) + idx = findfirst(==(net), c.nets) + @assert !isnothing(idx) + deleteat!(c.nets, idx) for (name, named_net) in c.net_names if named_net == net c.net_names[name] = nets[1] @@ -217,14 +219,14 @@ function topomat!(incidence::SparseMatrixCSC{T}) where {T<:Integer} row = 1; for col = 1:size(incidence)[2] - rows = filter(r -> r ≥ row, findall(!iszero, incidence[:, col])) + rows = filter(≥(row), findall(!iszero, incidence[:, col])) @assert length(rows) ≤ 2 isempty(rows) && continue t[col] = true; if rows[1] ≠ row - incidence[[rows[1], row], :] = incidence[[row, rows[1]], :] + incidence[[rows[1]::Integer, row], :] = incidence[[row, rows[1]], :] end if length(rows) == 2 @assert incidence[row, col] + incidence[rows[2], col] == 0 @@ -247,7 +249,7 @@ function topomat!(incidence::SparseMatrixCSC{T}) where {T<:Integer} tv = spzeros(T, size(dl, 2), size(incidence, 2)) # findall here works around JuliaLang/julia#27013 (introdcued in 0.7.0-DEV.4983) tv[:, findall(t)] = -dl' - tv[:, findall((!).(t))] = SparseMatrixCSC{T}(I, size(dl, 2), size(dl, 2)) + tv[:, findall((!).(t))] = SparseMatrixCSC{T}(I, size(dl, 2), size(dl, 2))::SparseMatrixCSC{T} tv, ti end @@ -467,9 +469,9 @@ ports. end tv, ti = topomat!(incid) - S = SparseMatrixCSC{Rational{BigInt}}([Mᵥ Mᵢ Mₓ Mₓ´ Mq; - blockdiag(tv, ti) spzeros(nb(circ)+numports, 2nx(circ)+nq(circ))]) - ũ, M = gensolve(S, SparseMatrixCSC{Rational{BigInt}}([Mu u0; spzeros(nb(circ)+numports, nu(circ)+1)])) + S = SparseMatrixCSC{Rational{BigInt}}([[Mᵥ Mᵢ Mₓ Mₓ´ Mq]; + [blockdiag(tv, ti) spzeros(nb(circ)+numports, 2nx(circ)+nq(circ))]]) + ũ, M = gensolve(S, SparseMatrixCSC{Rational{BigInt}}([[Mu u0]; spzeros(nb(circ)+numports, nu(circ)+1)])) # [v' i' x' xd' q']' = ũ + My for arbitrary y # now drop all rows concerning only internal voltages and currents indices = vcat(consecranges([nb(circ), numports, nb(circ), numports+2nx(circ)+nq(circ)])[[2,4]]...) diff --git a/src/kdtree.jl b/src/kdtree.jl index 1980f8d6..55723dfc 100644 --- a/src/kdtree.jl +++ b/src/kdtree.jl @@ -88,7 +88,7 @@ mutable struct Alts{T} end Alts(p::Vector{T}) where {T} = - Alts([AltEntry(1, zeros(T, length(p)), zero(T))], typemax(T), 0, 1) + Alts([AltEntry(1, zeros(T::Type, length(p)), zero(T))], typemax(T), 0, 1) function init!(alts::Alts{T}, best_dist, best_pidx) where {T} alts.number_valid = 1 diff --git a/src/solvers.jl b/src/solvers.jl index 0200d47a..7fdf6b60 100644 --- a/src/solvers.jl +++ b/src/solvers.jl @@ -136,6 +136,8 @@ function Base.copy!(dest::LinearSolver, src::LinearSolver) copyto!(dest.ipiv, src.ipiv) end +abstract type NonlinearSolver end + """ SimpleSolver @@ -146,7 +148,7 @@ the last solution found (or another solution provided externally) using the available Jacobians. Due to the missing global convergence, the `SimpleSolver` is rarely useful as such. """ -mutable struct SimpleSolver{NLEQ<:ParametricNonLinEq} +mutable struct SimpleSolver{NLEQ<:ParametricNonLinEq} <: NonlinearSolver nleq::NLEQ z::Vector{Float64} linsolver::LinearSolver @@ -242,7 +244,7 @@ can be combined with the `SimpleSolver` as `HomotopySolver{SimpleSolver}` to obtain a useful Newton homtopy solver with generally good convergence properties. """ -mutable struct HomotopySolver{BaseSolver} +mutable struct HomotopySolver{BaseSolver<:NonlinearSolver} <: NonlinearSolver basesolver::BaseSolver start_p::Vector{Float64} pa::Vector{Float64} @@ -314,7 +316,7 @@ See [M. Holters, U. Zölzer, "A k-d Tree Based Solution Cache for the Non-linear Equation of Circuit Simulations"](http://www.eurasip.org/Proceedings/Eusipco/Eusipco2016/papers/1570255150.pdf) for a more detailed discussion. """ -mutable struct CachingSolver{BaseSolver} +mutable struct CachingSolver{BaseSolver<:NonlinearSolver} <: NonlinearSolver basesolver::BaseSolver ps_tree::KDTree{Vector{Float64}, Matrix{Float64}} zs::Matrix{Float64} @@ -402,7 +404,7 @@ set_extrapolation_origin(solver::CachingSolver, p, z) = get_extrapolation_jacobian(solver::CachingSolver) = get_extrapolation_jacobian(solver.basesolver) -function linearize(solver, p::AbstractVector{Float64}) +function linearize(solver::NonlinearSolver, p::AbstractVector{Float64}) z = solve(solver, p) set_extrapolation_origin(solver, p, z) if !hasconverged(solver) diff --git a/test/runtests.jl b/test/runtests.jl index 43e11f44..169f1447 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,17 +7,17 @@ using ACME using FFTW: rfft using ProgressMeter: Progress, next! using SparseArrays: sparse, spzeros -using Test: @test, @test_broken, @test_logs, @test_throws, @testset +using Test: @inferred, @test, @test_broken, @test_logs, @test_throws, @testset @testset "topomat" begin - tv, ti = ACME.topomat(sparse([1 -1 1; -1 1 -1])) + tv, ti = @inferred ACME.topomat(sparse([1 -1 1; -1 1 -1])) @test tv*ti'==spzeros(2,1) # Pathological cases for topomat: # two nodes, one loop branch (short-circuited) -> voltage==0, current arbitrary - @test ACME.topomat(spzeros(Int, 2, 1)) == (hcat([1]), spzeros(0, 1)) + @test @inferred(ACME.topomat(spzeros(Int, 2, 1))) == (hcat([1]), spzeros(0, 1)) # two nodes, one branch between them -> voltage arbitrary, current==0 - @test ACME.topomat(sparse([1,2], [1,1], [1,-1])) == (spzeros(0, 1), hcat([1])) + @test @inferred(ACME.topomat(sparse([1,2], [1,1], [1,-1]))) == (spzeros(0, 1), hcat([1])) end @testset "LinearSolver" begin @@ -211,7 +211,7 @@ end @testset "gensolve/rank_factorize" begin a = Rational{BigInt}[1 1 1; 1 1 2; 1 2 1; 1 2 2; 2 1 1; 2 1 2] b = Rational{BigInt}[1 2 3 4 5 6; 6 5 4 3 2 1; 1 0 1 0 1 0] - nullspace = ACME.gensolve(sparse(a'), spzeros(Rational{BigInt}, size(a, 2), 0))[2] + nullspace = @inferred(ACME.gensolve(sparse(a'), spzeros(Rational{BigInt}, size(a, 2), 0)))[2] @test nullspace'*a == spzeros(3, 3) c, f = ACME.rank_factorize(sparse(a * b)) @test c*f == a*b @@ -234,13 +234,12 @@ end zuin in (zeros(Rational{BigInt}, 3, 1), hcat(Rational{BigInt}[1; 2; -1])) dq_full = p * dq + fq * zxin eq_full = p * eq + fq * zuin - mats = Dict{Symbol,Array}(:a => a, :b => b, :c => c, :dy => dy, :ey => ey, - :fy => fy, :dq_full => dq_full, :eq_full => eq_full, :fq => fq) - mats[:dq_fulls]=Matrix[mats[:dq_full]] - mats[:eq_fulls]=Matrix[mats[:eq_full]] - mats[:fqprev_fulls]=Matrix[mats[:eq_full]] - mats[:fqs]=Matrix[mats[:fq]] - ACME.reduce_pdims!(mats) + mats = ( + a = a, b = b, c = c, dy = dy, ey = ey, fy = fy, + dq_full = dq_full, eq_full = eq_full, fq = fq, + dq_fulls = [dq_full], eq_fulls = [eq_full], fqprev_fulls = [eq_full], fqs = [fq] + ) + mats = ACME.reduce_pdims!(mats) @test size(mats[:pexps][1], 2) == 3 @test mats[:pexps][1] * mats[:dqs][1] == mats[:dq_fulls][1] @test mats[:pexps][1] * mats[:eqs][1] == mats[:eq_fulls][1]