text:
"""
    emdsyn(sysm::MDMModel; rdim, nullspace = true, simple = false, minimal = true,
                           emtest = false, normalize = false, fast = true,
                           sdeg, smarg, poles, HDesign, MDtol, MDGainTol, MDfreq,
                           tcond, offset, atol, atol1, atol2, atol3, rtol)
                           -> (Q::MDFilter, R::MDFilterIF, info)

Solve the exact model detection problem (EMDP) for a given multiple synthesis model
`sysm::MDMModel`. The computed stable and proper filter objects `Q` and `R` contain the
model detection filter, representing the solution of the EMDP, and its internal form, respectively.

The returned named tuple `info`, with the components `info.tcond`, `info.degs`, `info.MDperf`, and `info.HDesign`,
contains additional synthesis related information (see below).

For a multiple synthesis model `sysm` containing `N` stable component models,
the `i`-th component system `sysm.sys[i]` has a descriptor realization of the form
`sysm.sys[i] = (Ai-λEi,Bi,Ci,Di)` and the corresponding input-output form is

      yi = Gui(λ)*u + Gdi(λ)*di + Gwi(λ)*wi + Gvi(λ)*vi

with the Laplace- or Z-transformed plant outputs `yi`, control inputs `u`,
disturbance inputs `di`, noise inputs `wi`, and auxiliary inputs `vi`,
and with `Gui(λ)`, `Gdi(λ)`, `Gwi(λ)` and
`Gvi(λ)`, the corresponding transfer function matrices.

The model detection filter object `Q`, contains in `Q.sys` the resulting bank of `N` filters.
The `i`-th filter `Q.sys[i]` is in a standard state-space form and generates `r_i`,
the `i`-th component (scalar or vector) of the overall residual vector `r := [r_1; r_2; ...; r_N]`.
The corresponding input-output (implementation) form of the `i`-th filter is

            r_i = Qyi(λ)*y + Qui(λ)*u   ,

where `Qyi(λ)` and `Qui(λ)` are the transfer function matrices from the output and control inputs to the `i`-th residual component.
The dimensions of output and control inputs are contained in the integers
`Q.ny` and `Q.mu`, respectively.

The model detection filter internal form object `R`, contains `R.sys`, the resulting array of `N × N`
internal form of the filters.
The `(i,j)`-th component filter `R.sys[i,j]` is in a standard state-space form,
generates the residual signal `r_ij`, and corresponds to the
input-output form

       r_ij = Ruij(λ)*u + Rdij(λ)*dj + Rwij(λ)*wj + Rvij(λ)*vj ,

where

       | Ruij(λ) Rdij(λ) Rwij(λ) Raij(λ) | := |Qyi(λ) Qui(λ)]*| Guj(λ) Gdj(λ) Gwj(λ) Gaj(λ) |.
                                                              |   I    0      0      0      |

The solution of the EMDP ensures that for the `i`-th filter, `Ruii(λ) = 0`, `Rdii(λ) = 0`, and
`[Ruij(λ) Rdij(λ)]` ``\\neq`` `0` for `j` ``\\neq`` `i`.

Various user options can be specified via keyword arguments as follows:

If `minimal = true` (default), least order filter synthesis is performed to determine
each of the component filters
`Q.sys[i]` for `i = 1, ...,N` and `R.sys[i,j]` for `i, j = 1, ...,N`  while
with `minimal = false` full order synthesis is performed.

If `HDesign = H`, then `H` is an `N`-dimensional array of full row rank or empty design matrices `H = [H_1, ..., H_N]`,
where `H_i` is the design matrix employed for the synthesis of the `i`-th component filter (default: `HDesign = missing`)

`rdim = q` specifies the vector `q`, whose `i`-th component `q[i]` specifies
the number of residual outputs for the `i`-th component filter `Q.sys[i]`.
If `q` is a scalar, then a vector `rdim` with all components equal to `q` is assumed.
The default value of `q[i]` is chosen as follows: if `HDesign = missing` or `H_i` is empty then
`q[i] = 1`, if `minimal = true`, or `q[i]` is the number of the nullspace basis
vectors used for the synthesis of `Q.sys [i]`, if `minimal = false`;
if  `H_i` specifies a full row rank design matrix, then `q[i]` is the row dimension of `H_i`.

If `emdtest = false` (default), only the input channels are used for model detectability tests.
If `emdtest = true`, extended model detectability tests are performed using both control and
disturbance input channels.

`MDfreq = ω` specifies a vector of real frequency values or a scalar real frequency value
for strong model detectability checks (default: `MDfreq = missing`).

If `nullspace = false` (default),  a full-order observer based nullspace basis is
used for the synthesis of the `i`-th filter whenever the `i`-th component system has no disturbance inputs.
Otherwise, minimal proper nullspace bases are used.
If `nullspace = true`, minimal proper nullspace bases are used for the synthesis of the
model detection filters.

If `simple = true`, simple proper nullspace bases are emplyed for synthesis.
The orders of the basis vectors employed for the synthesis of the `i`-th filter
are provided in `info.deg[i]`.
If `simple = false` (default), then no simple bases are computed.

`offset = β` specifies the boundary offset `β` to assess the stability of poles.
Accordingly, for the stability of a continuous-time system all real parts of poles must be at most `-β`,
while for the stability of a discrete-time system all moduli of poles must be at most `1-β`.
The default value used for `β` is `sqrt(ϵ)`, where `ϵ` is the working machine precision.

`smarg = α` specifies the stability margin which defines the stability
domain `Cs` of poles, as follows:
for a continuous-time system, `Cs` is the set of complex numbers
with real parts at most `α`,
while for a discrete-time system, `Cs` is the set of complex numbers with
moduli at most `α < 1` (i.e., the interior of a disc of radius `α` centered in the origin).
If `smarg` is missing, then the employed default values are `α = -β`
for a continuous-time system and `α = 1-β` for a discrete-time system,
where `β` is the boundary offset specified by the keyword argument `offset = β`.

`sdeg = γ` is the prescribed stability degree for the poles of the filters `Q` and `R`
(default: `γ = -0.05` for the real parts of poles for a continuous-time system and
`γ = 0.95` for the magnitudes of poles for a discrete-time system).

`poles = v` specifies a complex vector `v` containing a complex conjugate set
of desired poles within the stability domain `Cs` to be assigned for the filters `Q` and `R`
(default: `poles = missing`).

`tcond = tcmax` specifies the maximum alowed condition number `tcmax`
of the employed non-orthogonal transformations (default: `tcmax = 1.e4`).

`MDtol = tol1` specifies the threshold `tol1` for model detectability checks
   (default: `tol1 = 0.0001`).

`MDGainTol = tol2` specifies the threshold `tol2` for strong model detectability checks
   (default: `tol2 = 0.01`).

If `normalize = false` (default), the `i`-th component filter `Q.sys[i]` is scaled such that
the minimum gain of `R.sys[i,j]` for `j = 1, ..., N`,   `j` ``\\neq`` `i`, is equal to one.
If `normalize = true`, the standard normalization of component filters is performed to ensure
equal gains for `R.sys[1,j]` and `R.sys[j,1]`.

The rank determinations in the performed reductions
are based on rank revealing QR-decompositions with column pivoting
if `fast = true` (default) or the more reliable SVD-decompositions if `fast = false`.

The keyword arguments `atol1`, `atol2` and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of the state-space matrices
`Ai`, `Bi`, `Ci` and `Di`, the absolute tolerance for the nonzero elements of `Ei`,
and the relative tolerance for the nonzero elements of all above matrices.
The default relative tolerance is `ni*ϵ`, where `ϵ` is the working machine epsilon
and `ni` is the orders of the system `sysm.sys[i]`.
The keyword argument `atol3` is an absolute tolerance for observability tests
   (default: internally determined value).
The keyword argument `atol` can be used
to simultaneously set `atol1 = atol`, `atol2 = atol` and `atol3 = atol`.

The resulting named tuple `info` contains `(tcond, degs, MDperf, HDesign) `, where:

`info.tcond` is the maximum of the condition numbers of the employed
   non-orthogonal transformation matrices; a warning is issued if `info.tcond >= tcmax`;

`info.degs` is an `N`-dimensional vector, whose `i`-th component is an integer vector
containing the degrees of the basis vectors of the employed simple
nullspace basis for the synthesis of the `i`-th filter component, if `simple = true`,
and the degrees of the
basis vectors of an equivalent polynomial nullspace basis, if `simple = false`;

`info.MDperf` is an `N × N` array containing the resulting distance-mapping performance,
representing  the peak gains of the associated internal representations of
the model detection filters. If the `(i,j)`-th component filter `R.sys[i,j]`
has the input-output form

     rij = Ruij(λ)*u + Rdij(λ)*dj + Rwij(λ)*wj + Rvij(λ)*vj ,

then, the `(i,j)`-th performance gain is evaluated as
`info.MDperf[i,j] = ||Ruij(λ)||∞` if `emdtest = false` (default) or
`info.MDperf[i,j] = ||[Ruij(λ) Rdij(λ)]||∞` if `emdtest = true`.
If `MDfreq = ω`, then each gain `info.MDperf[i,j]` represents
the maximum of 2-norm pointwise gains evaluated for all frequencies in `ω`.

`info.HDesign` is an `N`-dimensional vector, whose `i`-th component
is the design matrix `H_i` employed for the synthesis of
the `i`-th model detection filter.

_Method:_ The Procedure EMD from [1] is implemented to solve
the exact model detection problem. For more details on
the least order synthesis of model detection filters see [2].

_References:_

[1] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques.
              Springer Verlag, 2017; sec. 5.2.

[2] A. Varga, Least order fault and model detection using multi-models.
       IEEE CDC'09, Shanghai, China, 2009.
"""
function emdsyn(sysm::MDMModel; rdim::Union{Vector{Int},Int,Missing} = missing, poles::Union{AbstractVector,Missing} = missing,
                      sdeg::Union{Real,Missing} = missing, smarg::Union{Real,Missing} = missing,
                      nullspace::Bool = false, minimal::Bool = true, simple::Bool = false, emdtest::Bool = false,
                      MDtol::Real = 0.0001, MDGainTol::Real = 0.01, MDfreq::Union{AbstractVector{<:Real},Real,Missing} = missing,
                      tcond::Real = 1.e4, HDesign::Union{Vector{<:Matrix},Matrix,Missing} = missing, normalize::Bool = false,
                      offset::Real = sqrt(eps(Float64)), atol::Real = zero(Float64), atol1::Real = atol, atol2::Real = atol, atol3::Real = atol,
                      rtol::Real = (size(sysm.sys[1].A,1)+1)*eps(Float64)*iszero(max(atol1,atol2)),
                      fast::Bool = true)
   N = length(sysm.sys)
   T1 = Float64
   # check stability of input model
   for i = 1:N
       isstable(sysm.sys[i]) || error("Input multiple model must be stable")
       #T1 = promote_type(T1,eltype(sysm.sys[i]))
   end
   p = size(sysm.sys[1],1);       # number of measurable outputs
   Ts = sysm.sys[1].Ts
   disc = (Ts != 0);  # system type (continuous- or discrete-time)

   # decode options

   strongMD = !ismissing(MDfreq)
   strongMD && !isa(MDfreq,Vector) && (MDfreq = [MDfreq])
   lfreq = strongMD ? length(MDfreq) : 0

   # set default stability degree
   sdegdefault = disc ? 0.95 : -0.05

   # stability margin
   ismissing(smarg) && (smarg = disc ? 1-offset : -offset)  # set default stability margin

   poles_nomissing = !ismissing(poles)

   # sort desired poles
   if poles_nomissing
      tempc = poles[imag.(poles) .> 0]
      if !isempty(tempc)
         tempc1 = conj(poles[imag.(poles) .< 0])
         isequal(tempc[sortperm(real(tempc))],tempc1[sortperm(real(tempc1))]) ||
                 error("poles must be a self-conjugated complex vector")
      end
      # check that all eigenvalues are inside of the stability region
      ( ((disc && any(abs.(poles) .> 1-offset) )  || (!disc && any(real.(poles) .> -offset)))  &&
            error("The elements of poles must lie in the stability region of interest") )
   end

   ismissing(sdeg) && ismissing(poles) && (sdeg = sdegdefault)  # set desired stability degree to default value

   # desired number of filter outputs
   if !ismissing(rdim)
      if isa(rdim,Vector)
         length(rdim) == N || error("dimension of rdim must be equal to the number of models")
         minimum(rdim) > 0 || error("all components of rdim must be positive")
      else
         rdim > 0 || error("rdim must be positive")
         rdim = fill(rdim, N)
      end
   end

   # imposed design option to form linear combinations of basis vectors
   emptyHD = ismissing(HDesign)

   if !emptyHD
       if isa(HDesign,Matrix)
          size(HDesign,1) == rank(HDesign) || error("HDesign must have full row rank")
          HDesign = fill(HDesign,N)
       else
          size(HDesign,1) == N || error("number of HDesign components must be equal to the number of models $N")
       end
       if !ismissing(rdim)
          for i = 1:N
              if !isempty(HDesign[i])
                 mH = size(HDesign[i],1)
                 if mH > 0
                    mH == rdim[i] || error("row dimension of HDesign[$i] must be equal to rdim[$i]")
                    mH == rank(HDesign[i]) || error("HDesign[$i] must have full row rank")
                 end
              end
          end
       end
   else
      HDesign = fill(T1[],N)
   end

   # decode input information
   mu = sysm.mu
   inpu = 1:mu
   md = sysm.md
   mw = sysm.mw
   ma = sysm.ma

   m = mu .+ (md+mw+ma)      # total number of inputs

   infotcond = similar(Vector{T1},N)
   infodegs = similar(Vector{Vector{Int}},N)
   HDesign1 = similar(Vector{Array{T1,2}},N)
   Qt = similar(Vector{DescriptorStateSpace{T1}},N)
   Rt = similar(Matrix{DescriptorStateSpace{T1}},N,N)
   distinf = fill(-1.,N,N)
   ismissing(rdim) ? rdimtarget = missing : rdimtarget = copy(rdim)
   for i = 1:N
      # solve the exact MD problem for the i-th system (with noise inputs)
      #
      # compute a minimal left nullspace basis Q1i for the i-th system
      desc = (sysm.sys[i].E != I)
      m2 = mw[i]+ma[i]
      sdegNS = strongMD ? sdegdefault : missing
      if nullspace || simple || md[i] > 0 || (desc && rcond(sysm.sys[i].E) < 1.e-7 )
         # form [ Gu Gd Gf Gw Gaux; I 0 0 0 0]
         syse = [sysm.sys[i]; eye(mu,m[i])];
         #
         # compute a left nullspace basis Q = Q1 of G1 = [Gu Gd; I 0] = 0 and
         # obtain QR = [ Q R ], where R = [ Rf Rw Raux] = Q*[Gf Gw Ga;0 0 0]
         qtemp, info1 = glnull(syse, m2; simple, atol1, atol2, rtol, fast, sdeg = sdegNS, offset)
         degs = info1.degs
         tcond1 = info1.tcond
      elseif mu == 0 && md[i] == 0
         qtemp = dss(eye(T1,p))
         degs = Int[]; tcond1 = 1.
      else
         # compute minimal basis as Q = Q1 = [ I -Gu] and set
         # QR = [ Q R ], where R = [ Gf Gw Ga ]
         qtemp = [ eye(p) dss(sysm.sys[i].A, sysm.sys[i].E, -sysm.sys[i].B[:,inpu],
                   sysm.sys[i].C, -sysm.sys[i].D[:,inpu]; Ts)]
         # perform stabilization if strong detectability has to be enforced
         degs = Int[]; tcond1 = 1.
      end
      infodegs[i] = degs

      nvec = size(qtemp,1);   # number of basis vectors
      # check solvability conditions
      nvec == 0 && error("emdsyn: empty nullspace basis - the EMDP is not solvable")

      nq = order(qtemp)             # order of the minimal basis

      # set H[i] for checking the solvability condition
      emptyHDi = isempty(HDesign[i]);

      if emptyHDi
         S = trues(nvec,max(1,lfreq),N)
         Htemp = eye(nvec)
         noutmax = nvec
      else
         mh, nh = size(HDesign[i])
         S = trues(mh,max(1,lfreq),N)
         degs = Int[];
         noutmax = mh
         if nh < nvec
            # pad with zeros: row rank is preserved
            Htemp = [ HDesign[i] zeros(T, mh, nvec-nh) ]
         else
            # remove trailing columns if necessary: rank may drop
            Htemp = HDesign[i][:,1:nvec];
            if nh > nvec && mh != rank(Htemp)
               error("The leading part of HDesign[$i] must have full row rank")
            end
         end
      end
      #strongMD && (qtemp = glcf(qtemp; atol1, atol2, rtol, fast, sdeg = sdegNS)[1])
      defer = !emptyHDi && minimal
      for j = 1:N
         # check j-th model detectability
         if i != j
            temp = gir(Htemp*qtemp*[sysm.sys[j]; eye(T1,mu,m[j])]; atol)
            #i == 1 && (println(size(temp)))
            if strongMD
               if emdtest
                  St = fdisspec_(temp[:,1:mu+md[j]], MDfreq; FDGainTol = MDGainTol, atol1, atol2, rtol = 0, fast, stabilize = false)[1]
               else
                  St = fdisspec_(temp[:,1:mu], MDfreq; FDGainTol = MDGainTol, atol1, atol2, rtol = 0, fast, stabilize = false)[1]
               end
               #i == 1 && (println("i= $i j = $j St = "); display(St))
               Sj = trues(size(St,1),lfreq)
               for k = 1:lfreq
                   Sj[:,k] = any(view(St,:,:,k),dims=2)
               end
               #i == 1 && (println("j = $j Sj = "); display(Sj))
               #println("i = $i j = $j Sj = "); display(Sj)
               if any(maximum(Sj,dims=1))
                  S[:,:,j] = Sj
               else
                  if defer
                     @warn "Checking of strong solvability condition #$i for model #$j deferred"
                  else
                    emdtest ? error("Strong extended model detection not feasible for model #$j") :
                              error("Strong model detection not feasible for model #$j ")
                  end
               end
            else
               if emdtest
                  Sj = fditspec_(temp[:,1:mu+md[j]]; atol1, atol2, rtol, FDtol = MDtol)
               else
                  Sj = fditspec_(temp[:,1:mu]; atol1, atol2, rtol, FDtol = MDtol)
               end
               if any(Sj)
                  S[:,:,j] = maximum(Sj,dims=2)
               else
                  if defer
                     @warn "Checking of solvability condition #$i for model #$j deferred"
                  else
                    emdtest ? error("Extended model detection not feasible for model #$j") :
                              error("Model detection not feasible for model #$j")
                  end
               end
            end
         end
      end
      #i == 1 && println("i = $i S = ")
      #display(S)
      #i == 1 && println("rdimtarget = $rdimtarget")
      # setup the number of filter outputs
      if ismissing(rdimtarget)
         if emptyHDi
            rdim = minimal ? 1 : nvec
         else
            rdim = size(HDesign[i],1)
         end
      else
         rdim = min(rdimtarget[i],nvec)
      end
      #i == 1 && println("rdim = $rdim")

      # Step 2): compute admissible Q2 to reduce the order of Q2*Q;
      # update Q <- Q2*Q, R <- Q2*R

      # reorder degs to correspond to the expected orders of basis vectors
      # corresponding to the actual order of outputs of QR
      reverse!(degs)
      if rdim < nvec
         #if rdim < noutmax && mf > 0
            # determine possible low order syntheses using i >= rmin basis vectors
            # and the corresponding expected orders

            finish = false    # set termination flag
            nout = rdim       # initialize number of selected basis vectors
            # if ~simple && minimal
            #    QR = xperm(QR,nq:-1:1);  # permute states to speedup glmcover1
            # end
            while !finish
                # choose nout basis vectors, which potentially lead to a least order
                # filter with rdim outputs:
                # basesel(i,:) contains the indices of candidate basis vectors;
                # ordsel(i)    contains the presumably achievable least orders
                #println("i = $i degs = $degs rdim = $rdim nout = $nout simple = $simple S = ")
                #i == 1 && display(S)
                basesel, ordsel = emdbasesel(S, degs, rdim, nout, simple)
                #
                # update the synthesis using the selections of candidate vector(s),
                # starting with the least (potentially) achievable order
                for ii = 1:size(basesel,1)
                    baseind = basesel[ii] # indices of current basis selection
                    if rdim == nout
                        hbase = eye(rdim)
                    else
                        hbase = rand(rdim,nout)
                    end
                    ip = [baseind; setdiff(1:nvec,baseind)][:]
                    if simple && !isempty(degs)
                       if minimal
                          # determine Q2i such that Q2i*Q1i is admisible and
                          # has least order
                          if emptyHDi
                             # select vectors and elliminate unobservable dynamics
                             noelim = falses(nq)
                             #ell = sum(degs(1:basesel(i,1)-1));
                             ell = sum(degs[1:basesel[ii][1]-1]);
                             for jj = 1:nout
                                 ellnext = sum(degs[1:baseind[jj]]);
                                 noelim[ell+1:ellnext] .= true;
                                 ell = ellnext
                             end
                          end
                          if rdim == nout
                             if emptyHDi
                                ir = noelim
                                Ar, Er, Br, Cr, Dr = dssdata(qtemp[baseind,:])
                                qtest = dss(view(Ar,ir,ir), Er == I ? I : view(Er,ir,ir),view(Br,ir,:),view(Cr,:,ir), Dr; Ts)
                                h = Htemp[ip[1:rdim],:];
                             else
                                qtest = gir(Htemp*qtemp; atol1, atol2, rtol);
                             end
                          else
                             # this case is possible only if HDesign[i] is empty
                             # build rdim linear combinations of the first nout vectors
                             ir = noelim
                             Ar, Er, Br, Cr, Dr = dssdata(qtest[baseind,:])
                             qtest = hbase*dss(view(Ar,ir,ir),Er == I ? I : view(Er,ir,ir),view(Br,ir,:),view(Cr,:,ir), Dr; Ts)
                             h = [ hbase zeros(T1, rdim, nvec-nout) ]
                             h = h[:,ip]  # permute columns to match unpermuted QR
                          end
                       else
                          # determine Q2i = Hi such that Q2i*Q1i is admisible
                          if rdim == nout
                             if emptyHDi
                                h = Htemp[ip[1:rdim],:]
                                qtest = gir(qtemp[baseind,:]; atol1, atol2, rtol, infinite = false)
                             else
                                qtest = gir(Htemp*qtemp; atol1, atol2, rtol, infinite = false)
                             end
                          else
                             # this case is possible only if HDesign[i] is empty
                             # build rdim linear combinations of the first nout vectors
                             qtest = gir(hbase*qtemp[baseind,:]; atol1, atol2, rtol, infinite = false)
                             h = [ hbase zeros(T1, rdim, nvec-nout) ]
                             h = h[:,ip]  # permute columns to match unpermuted QR
                          end
                       end
                    else
                       if minimal
                          # determine Q2i such that Q2i*Q1i is admisible and
                          # has least order
                          if rdim == nout
                             if emptyHDi
                                qtest, _, info2 = glmcover1(qtemp[ip,:], rdim; atol1, atol2, rtol)
                                if !isempty(ordsel) && (order(qtest) != ordsel[ii])
                                   @warn "emdsyn: expected reduced order not achieved"
                                end
                                h = Htemp[ip[1:rdim],:]
                                #i == 1 && println("h = $h")
                             else
                                qtest, _, info2 = glmcover1([Htemp; eye(nvec)]*qtemp[ip,:], rdim; atol1, atol2, rtol)
                             end
                          else
                             # this case is possible only if HDesign is empty
                             # build rdim linear combinations of the first nout vectors
                             h = [ hbase zeros(T1, rdim, nvec-nout) ];
                             qtest, _, info2 = glmcover1([h; eye(nvec)]*qtemp[ip,:], rdim; atol1, atol2, rtol)
                             h = h[:,ip]  # permute columns to match unpermuted QR
                           end
                       else
                          # determine Q2i = Hi such that Q2i*Q1i is admisible
                          if rdim == nout
                             if emptyHDi
                                h = Htemp[ip[1:rdim],:]
                                qtest = gir(qtemp[baseind,:]; atol1, atol2, rtol, infinite = false)
                             else
                                qtest = gir(Htemp*qtemp; atol1, atol2, rtol, infinite = false)
                             end
                          else
                             qtest = gir(hbase*qtemp[baseind,:]; atol1, atol2, rtol, infinite = false)
                             h = [ hbase zeros(T1, rdim, nvec-nout) ];
                             h = h[:,ip]  # permute columns to match unpermuted QR
                          end
                       end
                    end
                    # check model detectability of the current design;
                    if (rdim == nout && minimal) || rdim < nout
                       notOK = false
                       for j = 1:N
                           # check j-th model detectability
                           if i != j
                              # dismiss design if check fails
                              temp = gir(qtest*[sysm.sys[j]; eye(T1,mu,m[j])]; atol1, atol2, rtol)
                              if strongMD
                                 if emdtest
                                    St = fdisspec_(temp[:,1:mu+md[j]], MDfreq; FDGainTol = MDGainTol, atol1, atol2, rtol = 0, fast)[1]
                                 else
                                    St = fdisspec_(temp[:,1:mu], MDfreq; FDGainTol = MDGainTol, atol1, atol2, rtol = 0, fast)[1]
                                 end
                                 Sj = trues(size(St,1),lfreq);
                                 for k = 1:lfreq
                                     Sj[:,k] = any(view(St,:,:,k),dims=2)
                                 end
                                 notOK = !any(maximum(Sj,dims=1))
                                 #i == 1 && println("Sj = $Sj")
                                 #notOK = lfreq == 1 ? !any(St) : !all([any(view(St,:,:,k)) for k = 1:lfreq])
                              else
                                 if emdtest
                                    Sj = fditspec_(temp[:,1:mu+md[j]]; atol1, atol2, rtol, FDtol = MDtol)
                                 else
                                    Sj = fditspec_(temp[:,1:mu]; atol1, atol2, rtol, FDtol = MDtol)
                                 end
                                 #println("j = $j Sj = $Sj")
                                 notOK = !any(Sj)
                              end
                           end
                           notOK && break
                       end
                       if !notOK
                          if !simple && minimal
                             # adjust condition number of employed transformations
                             tcond1 = max(tcond1, info2.fnorm, info2.tcond)
                             tcond1 > tcond &&
                                @warn "emdsyn: possible loss of numerical stability due to ill-conditioned transformations"
         #                     if ~emptyHD
         #                        info.HDesign = Htemp;
         #                     end
                          end
                          qtemp = qtest
                          finish = true
                          break
                       end
                    else
                       qtemp = qtest
                       finish = true
                       break
                    end
                end
                if !finish
                   if emptyHDi
                      nout += 1
                      if nout > noutmax
                         finish = true
                         @warn "model detectability not achieved with the chosen number of residuals"
                      end
                   else
                     qtemp = gir(Htemp*qtemp; atol1, atol2, rtol);
                     minimal &&
                        (@warn "model detection with least order not feasible for the given tolerances")
                     finish = true;
                   end
                end
            end
            # set Q[i] = Q2i*Q1i
            Qt[i] = qtemp
            if emptyHDi
               Htemp = h;
            end
         else
            hbase = eye(T1,rdim)
            baseind = 1:nvec;
            # set Q[i] = Q2i*Q1i
            h = eye(T1,rdim)
            if emptyHDi
               Qt[i] = qtemp
               Htemp = h
            else
               Qt[i] = Htemp*qtemp
            end
         end

         # compute Q3i such that Q3i*Q[i] has a desired stability degree;
         # update Q[i] <- Q3i*Qi[i]
         k = 1;
         if simple && isequal(hbase,I) && emptyHDi
            # exploit the block diagonal structure of basis matrices al and cl
            # to compute block-diagonal Q3i
            al, el, bl, cl, dl, = dssdata(Qt[i])
            for ii = 1:length(baseind)
                blkord = degs[baseind[ii]]
                if blkord
                   i1 = k:k+blkord-1;
                   Qi = glcf(dss(al[i1,i1],el[i1,i1],bl[i1,:],cl[ii:ii,i1],dl[ii:ii,:];Ts);
                                atol1, atol2, atol3, sdeg, smarg, evals = poles)[1]
                   al[i1,i1] = Qi.A; bl[i1,:] = Qi.B;  cl[ii,i1] = Qi.C; dl[ii,:] = Qi.D
                   el[i1,i1] = (Qi.e == I) ? eye(T1,blkord) : Qi.E
                   k += blkord
                end
            end
            Qt[i] = dss(al, el, bl, cl, dl; Ts)
         else
            Qt[i] = glcf(Qt[i]; atol1, atol2, atol3, sdeg, smarg, evals = poles)[1]
         end

         for j = 1:N
           temp = gir(Qt[i]*[sysm.sys[j]; eye(T1,mu,m[j])]; atol1, atol2, rtol)
           # compute gains
           if i != j
              if strongMD
                 gains = zero(T1)
                 if emdtest
                    for k = 1:lfreq
                        gains = max(gains,norm(evalfr(temp[:,1:mu+md[j]]; fval = MDfreq[k])))
                    end
                 else
                    for k = 1:lfreq
                        gains = max(gains,opnorm(evalfr(temp[:,1:mu]; fval = MDfreq[k])))
                    end
                 end
                 distinf[i,j] = gains
              else
                 if emdtest
                    distinf[i,j] = glinfnorm(temp[:,1:mu+md[j]])[1]
                 else
                    distinf[i,j] = glinfnorm(temp[:,1:mu])[1]
                 end
              end
           else
              distinf[i,i] = 0
           end
           Rt[i,j] = temp
         end

         if normalize
            # scale to unit minimum gains
            scale = 1/minimum(view(distinf,i,[1:i-1; i+1:N]))
         else
            # scale to ensure distinf(i,1) = distinf(1,i)
            if i == 1
               scale = 1/minimum(view(distinf,1,2:N))
            else
               scale = distinf[1,i]/distinf[i,1]
            end
         end
         lmul!(scale,view(distinf,i,:))

         # scale and transform to standard state-space
         Qt[i]  = scale*gss2ss(Qt[i]; atol1, atol2, rtol)[1]
         for j = 1:N
            Rt[i,j] = scale*gss2ss(Rt[i,j]; atol1, atol2, rtol)[1]
         end
         infotcond[i] = tcond1
         HDesign1[i] = emptyHDi ? convert(Matrix{T1},h) : convert(Matrix{T1},Htemp)
   end

   Q = MDFilter{T1}(Qt, p, mu)
   #R = fdRset(QR[:,p+mu+1:end],faults = Vector(1:mf), noise = mf .+ Vector(1:mw), aux = (mf+mw) .+ Vector(1:ma))
   R = MDFilterIF{T1}(Rt, mu, md, mw, ma)
   info = (tcond = infotcond, degs = infodegs, HDesign = HDesign1, MDperf = distinf)

   return Q, R, info

   # end EMDSYN
end
function emdbasesel(S::BitArray, degs::Vector{Int}, rdim::Int, nout::Int, simple::Bool)
    #   emdbasesel(S, degs, rdim, nout, simple) -> (seli, selord)
    #
    #   Select admissible basis vectors for solving the exact model detection problem (EMDP)
    #   using the binary structure matrix `S` corresponding to `nvec` basis vectors.
    #   `S` is a `nvec x lfreq x N` binary 3-dimensional array or a `nvec x mf`
    #   binary 2-dimensional array (reshaped internally as an `nvec x mf x 1`
    #   binary 3-dimensional array).
    #   `seli` contains `nout`-touples of indices of basis vectors whose linear combination
    #   is admissible, i.e. , `S[seli[i],k,j]` has all elements nonzero for all `k = 1:lfreq`
    #   and `j = 1:N`. This ensures that the EMDP is solvable by using
    #   model detection filters with `rdim <= nout` outputs.
    #   If the associated `nvec` degrees contained in `degs` are provided, then
    #   `selord[i]` is the corresponding tentatively achievable least filter order.
    #   If `simple = true`, a simple basis is assumed, in which case, `degs[i]` is
    #   also the order of the `i`-th basis vector. If `simple = false`, a minimum
    #   rational basis is assumed. `selord` is empty if `degs` is empty.

    #   Method: The selection approach is used in conjunction with the synthesis
    #   Procedure EMD described in [1].

    #   References:
    #   [1] Varga A.
    #   Solving Fault Diagnosis Problems - Linear Synthesis Techniques. Springer Verlag, 2017.

    ndim = ndims(S)
    if ndim == 3
       nvec, lfreq, N = size(S)  # numbers of vectors, numbers of frequencies, numbers of models
    else
       nvec, N = size(S)   # numbers of vectors, number of models
       lfreq = 1
       S = reshape(S,nvec,1,N)
    end
    nd = length(degs)
    nodegs = (nd == 0)

    nodegs || length(degs) == nvec || error("the dimension of degs must be equal to the number of rows of S")

    (rdim >=1 && rdim <= nvec) || error("ndim must have a positive value not exceeding $nvec")
    (nout >=rdim && nout <= nvec) || error("nout must have a value at least $rdim and at most $nvec")

    nvec == 1 && (return [1], nodegs ? Int[] : degs )

    # find rdim combinations of nout vectors which solve the EFDP
    seli = collect(combinations(Vector(1:nvec),nout))
    ni = length(seli)
    selord = nodegs ? Int[] : fill(-1,ni)
    nqmax = sum(degs)
    ii = trues(ni)
    for i = 1:ni
        indv = seli[i];
        # check admissibility
        #stest = all(view(S,indv,:,:))
        stest = true;
        for j = 1:N
            #any(view(S,indv,:,j)) || (stest = false; break)
            for k = 1:lfreq
               stest = stest & any(view(S,indv,k,j))
            end
        end
        if stest
           if !nodegs
              # estimate orders
              if simple || rdim == nout
                 # degree = the sums of degrees of selected vectors
                 selord[i] = sum(degs[indv])
              else
                 # degree = rdim times the maximum degree of selected vectors
                 selord[i] = min(nqmax,rdim*maximum(degs[indv]))
              end
           end
        else
           ii[i] = false
        end
   end

    seli = seli[ii]
    if !nodegs
       selord = selord[ii];
       # sort row combinations to ensure increasing tentative orders
       ii = sortperm(selord)
       seli = seli[ii]
       selord = selord[ii];
    end
    return seli, selord
    # end EMDBASESEL
 end

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
import MathOptInterface

const MOI = MathOptInterface

mutable struct VariableInfo
    lower_bound::Float64
    upper_bound::Float64
end

VariableInfo() = VariableInfo(-Inf, Inf)

mutable struct Optimizer <: MOI.AbstractOptimizer
    result
    tol::Float64
    np:: Int64
    workers::Vector{Int64}
    nlp_data::MOI.NLPBlockData
    variable_info::Vector{VariableInfo}
    sense::MOI.OptimizationSense
    objective::Union{MOI.SingleVariable, MOI.ScalarAffineFunction{Float64}, MOI.ScalarQuadraticFunction{Float64}, Nothing}
    debug::Bool
end

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; copy_names = false)
    return MOI.Utilities.default_copy_to(model, src, copy_names)
end

struct EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator end

MOI.eval_objective(::EmptyNLPEvaluator, x) = NaN

empty_nlp_data() = MOI.NLPBlockData([], EmptyNLPEvaluator(), false)

function Optimizer(;workers = 2, tol = 1e-6, np = 30, debug = false)

    if workers > 1
        worker_ids = Distributed.workers()
        workers_present = nprocs()
        if  workers_present < 3     # True if the Julia session has only one worker
            addprocs(3 - workers_present)
            worker_ids = Distributed.workers()
        end
        @eval @everywhere using CharibdeOptim
        @eval @everywhere using JuMP
        @eval @everywhere using ModelingToolkit
        @eval @everywhere using IntervalArithmetic
        return Optimizer(nothing, tol, np, worker_ids, empty_nlp_data(), [], MOI.FEASIBILITY_SENSE, nothing, debug)
    else
        return Optimizer(nothing, tol, np, [1], empty_nlp_data(), [], MOI.FEASIBILITY_SENSE, nothing, debug)
    end
end

MOI.get(::Optimizer, ::MOI.SolverName) = "CharibdeOptim"

MOI.supports(::Optimizer, ::MOI.NLPBlock) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.SingleVariable}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.Interval{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.LessThan{Float64}}) = true

MOI.get(model::Optimizer, ::MOI.NumberOfVariables) = length(model.variable_info)

function MOI.get(model::Optimizer, ::MOI.ListOfVariableIndices)
    return [MOI.VariableIndex(i) for i in 1:length(model.variable_info)]
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveSense,
                 sense::MOI.OptimizationSense)
    model.sense = sense
    return
end

function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    model.nlp_data = nlp_data
    return
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction{F}, func::F) where {F <: Union{MOI.SingleVariable, MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction}}
    model.objective = func
    return
end

function MOI.is_empty(model::Optimizer)
    return isempty(model.variable_info) && model.nlp_data.evaluator isa EmptyNLPEvaluator && model.sense == MOI.FEASIBILITY_SENSE
end

function MOI.empty!(model::Optimizer)
    model.result = nothing
    empty!(model.variable_info)
    model.nlp_data = empty_nlp_data()
    model.sense = MOI.FEASIBILITY_SENSE
    model.objective = nothing
end

function MOI.add_constraint(model::Optimizer, var::MOI.SingleVariable, bound::MOI.Interval)
    vi = var.variable
    model.variable_info[vi.value].lower_bound = bound.lower
    model.variable_info[vi.value].upper_bound = bound.upper
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval}(vi.value)
end

function MOI.add_constraint(model::Optimizer, var::MOI.SingleVariable, hi::MOI.LessThan{Float64})
    vi = var.variable
    model.variable_info[vi.value].upper_bound = hi.upper
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}(vi.value)
end

function MOI.add_constraint(model::Optimizer, var::MOI.SingleVariable, low::MOI.GreaterThan{Float64})
    vi = var.variable
    model.variable_info[vi.value].lower_bound = low.lower
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}(vi.value)
end

function MOI.add_variable(model::Optimizer)
    push!(model.variable_info, VariableInfo())
    return MOI.VariableIndex(length(model.variable_info))
end
function MOI.add_variables(model::Optimizer, n::Int)
    return [MOI.add_variable(model) for i in 1:n]
end

function eval_function(var::MOI.SingleVariable, x)
    return x[var.variable.value]
end

function eval_function(aff::MOI.ScalarAffineFunction, x)
    function_value = aff.constant
    for term in aff.terms
        function_value += term.coefficient*x[term.variable_index.value]
    end
    return function_value
end

function eval_function(quad::MOI.ScalarQuadraticFunction, x)
    function_value = quad.constant
    for term in quad.affine_terms
        function_value += term.coefficient*x[term.variable_index.value]
    end
    for term in quad.quadratic_terms
        row_idx = term.variable_index_1
        col_idx = term.variable_index_2
        coefficient = term.coefficient
        if row_idx == col_idx
            function_value += 0.5*coefficient*x[row_idx.value]*x[col_idx.value]
        else
            function_value += coefficient*x[row_idx.value]*x[col_idx.value]
        end
    end
    return function_value
end

function substitute_variables(expr::Expr)
    if expr.head == :ref && length(expr.args) == 2 && expr.args[1] == :x
        return :(x[$(expr.args[2].value)])
    else
        for (index, arg) in enumerate(expr.args)
            expr.args[index] = substitute_variables(arg)
        end
    end
    return expr
end

substitute_variables(arg) = arg

function substitute_symbols(expr::Expr)
    if expr.head == :ref && length(expr.args) == 2 && expr.args[1] == :x
        return Symbol(:(x[$(expr.args[2].value)]))
    else
        for (index, arg) in enumerate(expr.args)
            expr.args[index] = substitute_symbols(arg)
        end
    end
    return expr
end
substitute_symbols(arg) = arg

function preprocess!(ex)
    isa(ex, Symbol) && return Variable(ex; known=false)()
    if isa(ex, Expr) && ex.head === :call
        ex.args[1] = eval(ex.args[1])
        @. ex.args[2:end] = preprocess!(ex.args[2:end])
    end
    return ex
end

function eval_objective(model::Optimizer, eval_expr::Union{Function, Nothing}, x)
    if eval_expr !== nothing
        return invokelatest(eval_expr, x)
    elseif model.objective !== nothing
        return eval_function(model.objective, x)
    else
        return 0.0
    end
end

function parse_constraint(ex::Expr)
    expr =
    @match ex begin
       ((a_ <= b_) | (a_ < b_) | (a_ ≤ b_)) => (a)
       ((a_ >= b_) | (a_ > b_) | (a_ ≥ b_)) => (a)
       ((a_ == b_) | (a_ = b_))   => (a)
       ((a_ <= b_ <= c_) | (a_ < b_ < c_) | (a_ <= b_ < c) | (a_ < b_ <= c)) => (b)
       ((a_ >= b_ >= c_) | (a_ > b_ > c_) | (a_ >= b_ > c_) | (a_ > b_ >= c)) => (b)
       ((a_ ∈ [b_, c_]) | (a_ in [b_, c_]) | (a_ ∈ b_ .. c_) | (a_ in b_ .. c_))  => (a)
   end
   return (expr)
end

function construct_constraints(model::Optimizer, num_constraints)
    vars = [Variable(Symbol(:(x[$i])))() for i in 1:length(model.variable_info)]
    constraints = CharibdeOptim.Constraint{Float64}[]

    for i in 1:num_constraints
        expr = MOI.constraint_expr(model.nlp_data.evaluator, i)
        expr = parse_constraint(expr)
        expr = substitute_symbols(expr)
        operation = convert(ModelingToolkit.Expression, preprocess!(expr))
        push!(constraints, constraint(vars, operation, Interval(model.nlp_data.constraint_bounds[i].lower, model.nlp_data.constraint_bounds[i].upper)))
    end
    return constraints
end

function diffevol_worker(model::Optimizer, search_space::IntervalBox{N,T}, ch_master_to_slave::RemoteChannel{Channel{Tuple{IntervalBox{N, T}, T}}}, ch_slave_to_master::RemoteChannel{Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}} ) where{N,T}

    eval_expr = nothing
    if model.nlp_data.has_objective
        MOI.initialize(model.nlp_data.evaluator, [:ExprGraph])
        expr = MOI.objective_expr(model.nlp_data.evaluator)
        expr = substitute_variables(expr)
        eval_expr = eval(:(x -> $(expr)))
    end

    obj_func(x) = eval_objective(model, eval_expr, x)
    num_constraints = length(model.nlp_data.constraint_bounds)

    if num_constraints == 0
        if model.sense == MOI.MIN_SENSE
            diffevol_minimise(obj_func, search_space, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, np = model.np, de_ind = false)
        else
            @assert sense == MOI.MAX_SENSE
            diffevol_maximise(obj_func, search_space, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, np = model.np, de_ind = false)
        end
    else
        constraints = construct_constraints(model, num_constraints)
        if model.sense == MOI.MIN_SENSE
            diffevol_minimise(obj_func, search_space, constraints, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, np = model.np, de_ind = false)
        else
            @assert model.sense == MOI.MAX_SENSE
            diffevol_maximise(obj_func, search_space, constraints, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, np = model.np, de_ind = false)
        end
    end
end

function ibc_worker(model::Optimizer, search_space::IntervalBox{N,T}, debug::Bool, ch_master_to_slave::RemoteChannel{Channel{Tuple{IntervalBox{N, T}, T}}}, ch_slave_to_master::RemoteChannel{Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}}, tol::T) where{N,T}

    eval_expr = nothing
    if model.nlp_data.has_objective
        MOI.initialize(model.nlp_data.evaluator, [:ExprGraph])
        expr = MOI.objective_expr(model.nlp_data.evaluator)
        expr = substitute_variables(expr)
        eval_expr = eval(:(x -> $(expr)))
    end

    obj_func(x) = eval_objective(model, eval_expr, x)

    num_constraints = length(model.nlp_data.constraint_bounds)

    if num_constraints == 0
        if model.sense == MOI.MIN_SENSE
            ibc_minimise(obj_func, search_space, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, debug = debug, tol = tol, ibc_ind = false)
        else
            @assert sense == MOI.MAX_SENSE
            ibc_maximise(obj_func, search_space, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, debug = debug, tol = tol, ibc_ind = false)
        end
    else
        constraints = construct_constraints(model, num_constraints)
        if model.sense == MOI.MIN_SENSE
            ibc_minimise(obj_func, search_space, constraints, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, debug = debug, tol = tol, ibc_ind = false)
        else
            @assert model.sense == MOI.MAX_SENSE
            ibc_maximise(obj_func, search_space, constraints, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, debug = debug, tol = tol, ibc_ind = false)
        end
    end

end

function optimize_serial(model::Optimizer, obj_func::Function, search_space::IntervalBox{N,T}) where{N,T}
    ch_master_to_slave = Channel{Tuple{IntervalBox{N, T}, T}}(1)
    ch_slave_to_master = Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}(1)

    num_constraints = length(model.nlp_data.constraint_bounds)

    if num_constraints == 0
        if model.sense == MOI.MIN_SENSE
            r1 = @async diffevol_minimise(obj_func, search_space, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, np = model.np, de_ind = false)
            r2 = @async ibc_minimise(obj_func, search_space, debug = model.debug, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, tol = model.tol, ibc_ind = false)
            model.result = fetch(r2)
        else
            @assert model.sense == MOI.MAX_SENSE
            r1 = @async diffevol_maximise(obj_func, search_space, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, np = model.np, de_ind = false)
            r2 = @async ibc_maximise(obj_func, search_space, debug = model.debug, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, tol = model.tol, ibc_ind = false)
            model.result = fetch(r2)
        end
    else
        constraints = construct_constraints(model, num_constraints)
        if model.sense == MOI.MIN_SENSE
            r1 = @async diffevol_minimise(obj_func, search_space, constraints, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, np = model.np, de_ind = false)
            r2 = @async ibc_minimise(obj_func, search_space, constraints, debug = model.debug, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, tol = model.tol, ibc_ind = false)
            model.result = fetch(r2)

        else
            @assert model.sense == MOI.MAX_SENSE
            r1 = @async diffevol_maximise(obj_func, search_space, constraints, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, np = model.np, de_ind = false)
            r2 = @async ibc_maximise(obj_func, search_space, constraints, debug = model.debug, ibc_chnl = ch_master_to_slave, diffevol_chnl = ch_slave_to_master, tol = model.tol, ibc_ind = false)
            model.result = fetch(r2)
        end

    end
end

function optimize_parallel(model::Optimizer, obj_func::Function, search_space::IntervalBox{N,T}) where{N,T}

    ch_master_to_slave = RemoteChannel(()->Channel{Tuple{IntervalBox{N, T}, T}}(1))
    ch_slave_to_master = RemoteChannel(()->Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}(1))

    r1 = remotecall(CharibdeOptim.diffevol_worker, model.workers[1], model, search_space, ch_master_to_slave, ch_slave_to_master)
    r2 = remotecall(CharibdeOptim.ibc_worker, model.workers[2], model, search_space, model.debug, ch_master_to_slave, ch_slave_to_master, model.tol)
    model.result = fetch(r2)
end

function MOI.optimize!(model::Optimizer)

    eval_expr = nothing
    if model.nlp_data.has_objective
        MOI.initialize(model.nlp_data.evaluator, [:ExprGraph])
        expr = MOI.objective_expr(model.nlp_data.evaluator)
        expr = substitute_variables(expr)
        eval_expr = eval(:(x -> $(expr)))
    end

    obj_func(x) = eval_objective(model, eval_expr, x)

    X = [Interval(var.lower_bound, var.upper_bound) for var in model.variable_info]
    search_space = IntervalBox(X...)

    if model.workers[1] != 1
        optimize_parallel(model, obj_func, search_space)
    else
        optimize_serial(model, obj_func, search_space)
    end
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    if model.result == nothing
        return MOI.OPTIMIZE_NOT_CALLED
    else
        return MOI.OPTIMAL
    end
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveValue)
    if model.result === nothing
        error("ObjectiveValue not available.")
    end
    return model.result[1].lo
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return (model.result !== nothing) ? 1 : 0
end

function MOI.get(model::Optimizer, ::MOI.PrimalStatus)
    if model.result === nothing
        return MOI.NO_SOLUTION
    else
        return MOI.FEASIBLE_POINT
    end
end

function MOI.get(model::Optimizer, ::MOI.DualStatus)
        return MOI.NO_SOLUTION
end

function MOI.get(model::Optimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex)
    if model.result === nothing
        error("VariablePrimal not available.")
    end
    return model.result[2][1][vi.value].lo
end

ext: jl

max_issues_repo_name: yashcodes/Charibde

max_forks_repo_name: yashcodes/Charibde


text:
using Plots
using LinearAlgebra
using SparseArrays
using Arpack
using ForwardDiff
using BenchmarkTools
using Formatting
using FileIO

#readfile
println("opening file... \n")
a = open("../programing/efelkismos-simple.unv","r")
include("readfile.jl") #returns array of: nodes, edge, elem, bc

include("elastic_functions.jl") #defines basic functions

#start building global array
n = length(elem[:,1])
dofs = 2*length(nodes[:,1])
global E = Ematrix(220e9,0.5,"stress")
global K = sparse([],[],Float64[],dofs,dofs);
global M = sparse([],[],Float64[],dofs,dofs);

println("calculating stiffness/mass matrix...\n")
h =1
@time begin
for i =1:n
    i1 = elem[i,:]'
    xs = nodes[i1,1]
    ys = nodes[i1,2]
    i1 = [2 .*i1 (i1 .* 2 .-1)]
    i2 = sort(i1,dims=2)
    @show i2
    kre = h .* ke(E,xs,ys)
    mre = h .* 7e3 .* me(xs,ys)
    K[i1,i1] = K[i1,i1][1,:,1,:] + kre
    M[i1,i1] = M[i1,i1][1,:,1,:] + mre
end
end

include("rhs.jl") #calculate F

#gib traction nodes where (u,v) are known
#reminder, ix contains boundary condition nodes INDEXES
println("solving for U....\n")
freedofs = zeros(dofs,1)
fixedofs = edge[bc[ix[1]],:][:]'

freedofs[[2 .*fixedofs (2 .* fixedofs .-1)]] .= 1
freedofs = [i for i=1:dofs if freedofs[i]==0]
fixedofs = [2 .* fixedofs (2 .* fixedofs .-1)]
U = zeros(dofs,1)
rhs = F
rhs = rhs[freedofs] - K[freedofs,fixedofs'][:,:,1]*U[fixedofs']
Kfree = K[freedofs,freedofs]
Mfree = M[freedofs,freedofs]

Ufree = U[freedofs]
U[freedofs] = Kfree\rhs

n = length(elem[:,1])
for i=1:n
    i1 = elem[i,:]'
    xs = nodes[i1,1]
    ys = nodes[i1,2]
    i1 = [2 .*i1 (i1 .* 2 .-1)]
end

println("calculating element average stress....\n")
#calculate element stress
include("stress.jl")
#println(stress')
println("DONE! \n")

i1 = elem[1,:]'
xs = nodes[i1,1]
ys = nodes[i1,2]

println(E*b([1 1 1]./3,xs,ys)*U[[2 .* i1 (2 .* i1 .- 1)]]')

xs = [1 0 0 0.5 0 0.5]
ys = [0 -.5 .5 -.25 0 .25]
l = [1 0 0]

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
##### update physiological attributes of each individual
function plankton_growth!(plank, nuts, rnd, p, ΔT, t, arch::Architecture)
    calc_inorganic_uptake!(plank, nuts, p, arch)

    update_quotas_1!(plank, ΔT, arch)

    calc_organic_uptake!(plank, nuts, p, arch)

    calc_ρChl!(plank, nuts.par, p, arch)

    calc_respir!(plank, nuts.T, p, arch)

    update_quotas_2!(plank, ΔT, p, arch)

    calc_BS!(plank, p, arch)
    update_biomass!(plank, p, ΔT, arch)
    update_cellsize!(plank, p, arch)

    ##### probabilities of grazing, mortality, and cell division
    ##### check the probabilities every 10 time steps or 1 hour whichever is shorter
    if t%(ΔT*(min(10,3600÷ΔT))) == 0
        calc_graz_quadratic!(plank, nuts, p.grz_P, arch)
        calc_mort!(plank, p, arch)
        calc_dvid!(plank, divide_type(p.dvid_type), p, t, arch)
        get_probability!(plank, rnd, ΔT, arch)
    else
        @inbounds plank.graz .= 0.0
        @inbounds plank.mort .= 0.0
        @inbounds plank.dvid .= 0.0
    end
end


text:
N=2;
A=CudaArray(rand(N,N));
B=CudaArray(rand(N,N));
C=CudaArray(rand(N,N));

        device(0)

for i=1:10
tic()
#result = devices(dev->capability(dev)[1]>=2) do devlist
#    mymodule.init(devlist) do dev
        vadd(A,B,C)
        #to_host(C)-to_host(B)-to_host(A)
#    end
end
toc()
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
module LibFEM
"""
    function declaration: d1_spring_assemble(K,k,i,j)

This function assembles the element stiffness
matrix k of the spring with nodes i & j into the
global stiffness matrix K.
This function returns the global stiffness matrix K
after the element stiffness matrix k is assembled.
"""
function d1_spring_assemble(K, k, i, j)
    K[i, i] = K[i, i] + k[1, 1]
    K[i, j] = K[i, j] + k[1, 2]
    K[j, i] = K[j, i] + k[2, 1]
    K[j, j] = K[j, j] + k[2, 2]
    return K
end
export d1_spring_assemble
"""
    function declaration: d1_spring_elementforce(k,u)

This function returns the element nodal force
vector given the element stiffness matrix k
and the element nodal displacement vector u.
"""
    function d1_spring_elementforce(k, u)
    return k * u
end
export d1_spring_elementforce
"""
    function declaration: d1_spring_elementstiffness(k)

This function returns the element stiffness
matrix for a spring with stiffness k.
The size of the element stiffness matrix
is 2 x 2.
"""
function d1_spring_elementstiffness(k)
    return [k -k; -k k]
end
export d1_spring_elementstiffness
"""
    function declaration: d1_truss_assemble(K,k,i,j)

This function assembles the element stiffness
matrix k of the linear bar with nodes i & j
into the global stiffness matrix K.
This function returns the global stiffness
matrix K after the element stiffness matrix
k is assembled.
"""
function d1_truss_assemble(K, k, i, j)
    K[i, i] = K[i, i] + k[1, 1]
    K[i, j] = K[i, j] + k[1, 2]
    K[j, i] = K[j, i] + k[2, 1]
    K[j, j] = K[j, j] + k[2, 2]
    return K
end
export d1_truss_assemble
"""
    function declaration: d1_truss_elementforce(k,u)

This function returns the element nodal
force vector given the element stiffness
matrix k & the element nodal displacement
vector u.
"""
function d1_truss_elementforce(k, u)
    return k * u
end
export d1_truss_elementforce
"""
    function declaration: d1_truss_elementstiffness(E,A,L)

This function returns the element
stiffness matrix for a linear bar with
modulus of elasticity E; cross-sectional
area A; & length L. The size of the
element stiffness matrix is 2 x 2.
"""
function d1_truss_elementstiffness(E, A, L)
    return [E * A / L -E * A / L; -E * A / L E * A / L]
end
export d1_truss_elementstiffness
"""
    function declaration: d1_truss_elementstress(k,u,A)

This function returns the element nodal
stress vector given the element stiffness
matrix k; the element nodal displacement
vector u; & the cross-sectional area A.
"""
function d1_truss_elementstress(k, u, A)
    return k * u / A
end
export d1_truss_elementstress
"""
    function declaration: d1_truss_elementstrain(k,u,A)

This function returns the element nodal
stress vector given the element stiffness
matrix k; the element nodal displacement
vector u; & the cross-sectional area A.
"""
function d1_truss_elementstrain(L, u)
    return 1/L * u
end
export d1_truss_elementstrain
"""
    function declaration: d2_beam_assemble(K,k,i,j)

This function assembles the element stiffness
matrix k of the plane beam element with nodes
i & j into the global stiffness matrix K.
This function returns the global stiffness
matrix K after the element stiffness matrix
k is assembled.
"""
function d2_beam_assemble(K, k, i, j)
    K[3*i-2, 3*i-2] = K[3*i-2, 3*i-2] + k[1, 1]
    K[3*i-2, 3*i-1] = K[3*i-2, 3*i-1] + k[1, 2]
    K[3*i-2, 3*i] = K[3*i-2, 3*i] + k[1, 3]
    K[3*i-2, 3*j-2] = K[3*i-2, 3*j-2] + k[1, 4]
    K[3*i-2, 3*j-1] = K[3*i-2, 3*j-1] + k[1, 5]
    K[3*i-2, 3*j] = K[3*i-2, 3*j] + k[1, 6]
    K[3*i-1, 3*i-2] = K[3*i-1, 3*i-2] + k[2, 1]
    K[3*i-1, 3*i-1] = K[3*i-1, 3*i-1] + k[2, 2]
    K[3*i-1, 3*i] = K[3*i-1, 3*i] + k[2, 3]
    K[3*i-1, 3*j-2] = K[3*i-1, 3*j-2] + k[2, 4]
    K[3*i-1, 3*j-1] = K[3*i-1, 3*j-1] + k[2, 5]
    K[3*i-1, 3*j] = K[3*i-1, 3*j] + k[2, 6]
    K[3*i, 3*i-2] = K[3*i, 3*i-2] + k[3, 1]
    K[3*i, 3*i-1] = K[3*i, 3*i-1] + k[3, 2]
    K[3*i, 3*i] = K[3*i, 3*i] + k[3, 3]
    K[3*i, 3*j-2] = K[3*i, 3*j-2] + k[3, 4]
    K[3*i, 3*j-1] = K[3*i, 3*j-1] + k[3, 5]
    K[3*i, 3*j] = K[3*i, 3*j] + k[3, 6]
    K[3*j-2, 3*i-2] = K[3*j-2, 3*i-2] + k[4, 1]
    K[3*j-2, 3*i-1] = K[3*j-2, 3*i-1] + k[4, 2]
    K[3*j-2, 3*i] = K[3*j-2, 3*i] + k[4, 3]
    K[3*j-2, 3*j-2] = K[3*j-2, 3*j-2] + k[4, 4]
    K[3*j-2, 3*j-1] = K[3*j-2, 3*j-1] + k[4, 5]
    K[3*j-2, 3*j] = K[3*j-2, 3*j] + k[4, 6]
    K[3*j-1, 3*i-2] = K[3*j-1, 3*i-2] + k[5, 1]
    K[3*j-1, 3*i-1] = K[3*j-1, 3*i-1] + k[5, 2]
    K[3*j-1, 3*i] = K[3*j-1, 3*i] + k[5, 3]
    K[3*j-1, 3*j-2] = K[3*j-1, 3*j-2] + k[5, 4]
    K[3*j-1, 3*j-1] = K[3*j-1, 3*j-1] + k[5, 5]
    K[3*j-1, 3*j] = K[3*j-1, 3*j] + k[5, 6]
    K[3*j, 3*i-2] = K[3*j, 3*i-2] + k[6, 1]
    K[3*j, 3*i-1] = K[3*j, 3*i-1] + k[6, 2]
    K[3*j, 3*i] = K[3*j, 3*i] + k[6, 3]
    K[3*j, 3*j-2] = K[3*j, 3*j-2] + k[6, 4]
    K[3*j, 3*j-1] = K[3*j, 3*j-1] + k[6, 5]
    K[3*j, 3*j] = K[3*j, 3*j] + k[6, 6]
    return K
end
export d2_beam_assemble
"""
    function declaration: d2_beam_elementaxialdiagram(f, L)

This function plots the axial force
diagram for the plane beam element
with nodal force vector f & length L.
"""
function d2_beam_elementaxialdiagram(f, L)
    x = [0 L]'
    z = [-f[1] f[4]]'
    #hold on
    title("Axial Force Diagram")
    plot(x, z)
    y1 = [0 0]'
    plot(x, y1, 'k')
end
export d2_beam_elementaxialdiagram
"""
    function declaration d2_beam_elementforce(E,A,I,L,theta,u)

This function returns the element force
vector given the modulus of elasticity E
the cross-sectional area A; the moment of
inertia I; the length L; the angle theta
(in degrees), & the element nodal
displacement vector u.
"""
function d2_beam_elementforce(E, A, I, L, theta, u)
    x = theta * pi / 180
    C = cos(x)
    S = sin(x)
    w1 = E * A / L
    w2 = 12 * E * I / (L * L * L)
    w3 = 6 * E * I / (L * L)
    w4 = 4 * E * I / L
    w5 = 2 * E * I / L
    kprime = [
        w1 0 0 -w1 0 0
        0 w2 w3 0 -w2 w3
        0 w3 w4 0 -w3 w5
        -w1 0 0 w1 0 0
        0 -w2 -w3 0 w2 -w3
        0 w3 w5 0 -w3 w4
    ]
    T = [
        C S 0 0 0 0
        -S C 0 0 0 0
        0 0 1 0 0 0
        0 0 0 C S 0
        0 0 0 -S C 0
        0 0 0 0 0 1
    ]
    return kprime * T * u
end
export d2_beam_elementforce
"""
    function declaration: d2_beam_elementlength(x1,y1,x2,y2)

This function returns the length of the
plane beam element whose first node has
coordinates [x1,y1] & second node has
coordinates [x2,y2].
"""
function d2_beam_elementlength(x1, y1, x2, y2)
    return sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1))
end
export d2_beam_elementlength
"""
    function declaration: d2_beam_elementmomentdiagram(f, L)

This function plots the bending
moment diagram for the plane beam
element with nodal force vector f & length L.
"""
function d2_beam_elementmomentdiagram(f, L)
    x = [0 L]'
    z = [-f[3] f[6]]'
    #hold on
    title("Bending Moment Diagram")
    plot(x, z)
    y1 = [0 0]'
    plot(x, y1, 'k')
end
export d2_beam_elementmomentdiagram
"""
    function declaration: d2_beam_elementsheardiagram(f, L)

This function plots the shear force
diagram for the plane beam element
with nodal force vector f & length L.
"""
function d2_beam_elementsheardiagram(f, L)
    x = [0 L]'
    z = [f[2] -f[5]]'
    #hold on
    title("Shear Force Diagram")
    plot(x, z)
    y1 = [0 0]'
    plot(x, y1, 'k')
end
export d2_beam_elementsheardiagram
"""
    function declaration: d2_beam_elementstiffness(E,A,I,L,theta)

This function returns the element
stiffness matrix for a plane beam
element with modulus of elasticity E;
cross-sectional area A; moment of
inertia I; length L; & angle()
theta [in degrees].
The size of the element stiffness matrix is 6 x 6.
"""
function d2_beam_elementstiffness(E, A, I, L, theta)
    x = theta * pi / 180
    C = cos(x)
    S = sin(x)
    w1 = A * C * C + 12 * I * S * S / (L * L)
    w2 = A * S * S + 12 * I * C * C / (L * L)
    w3 = (A - 12 * I / (L * L)) * C * S
    w4 = 6 * I * S / L
    w5 = 6 * I * C / L
    return E / L * [
        w1 w3 -w4 -w1 -w3 -w4
        w3 w2 w5 -w3 -w2 w5
        -w4 w5 4 * I w4 -w5 2 * I
        -w1 -w3 w4 w1 w3 w4
        -w3 -w2 -w5 w3 w2 -w5
        -w4 w5 2 * I w4 -w5 4 * I
    ]
end
export d2_beam_elementstiffness
"""
    function declaration:  d2_spring_assemble(K,k,i,j)

This function assembles the element stiffness
matrix k of the 2D spring element with nodes
i & j into the global stiffness matrix K.
This function returns the global stiffness
matrix K after the element stiffness matrix
k is assembled.
"""
function d2_spring_assemble(K, k, i, j)
    K[2*i-1, 2*i-1] = K[2*i-1, 2*i-1] + k[1, 1]
    K[2*i-1, 2*i] = K[2*i-1, 2*i] + k[1, 2]
    K[2*i-1, 2*j-1] = K[2*i-1, 2*j-1] + k[1, 3]
    K[2*i-1, 2*j] = K[2*i-1, 2*j] + k[1, 4]
    K[2*i, 2*i-1] = K[2*i, 2*i-1] + k[2, 1]
    K[2*i, 2*i] = K[2*i, 2*i] + k[2, 2]
    K[2*i, 2*j-1] = K[2*i, 2*j-1] + k[2, 3]
    K[2*i, 2*j] = K[2*i, 2*j] + k[2, 4]
    K[2*j-1, 2*i-1] = K[2*j-1, 2*i-1] + k[3, 1]
    K[2*j-1, 2*i] = K[2*j-1, 2*i] + k[3, 2]
    K[2*j-1, 2*j-1] = K[2*j-1, 2*j-1] + k[3, 3]
    K[2*j-1, 2*j] = K[2*j-1, 2*j] + k[3, 4]
    K[2*j, 2*i-1] = K[2*j, 2*i-1] + k[4, 1]
    K[2*j, 2*i] = K[2*j, 2*i] + k[4, 2]
    K[2*j, 2*j-1] = K[2*j, 2*j-1] + k[4, 3]
    K[2*j, 2*j] = K[2*j, 2*j] + k[4, 4]
    return K
end
export  d2_spring_assemble
"""
    function declaration: d2_spring_elementforce(k,theta,u)

This function returns the element force
given the stiffness k &
the angle theta [in degrees], and the
element nodal displacement vector u.
"""
function d2_spring_elementforce(k, theta, u)
    x = theta * pi / 180
    C = cos(x)
    S = sin(x)
    return k * [-C -S C S] * u
end
export d2_spring_elementforce
"""
    function declaration: d2_spring_elementstiffness(k,theta)

This function returns the element
stiffness matrix for a 2D spring
with stiffness k &
angle theta [in degrees].
The size of the element stiffness
matrix is 4 x 4.
"""
function d2_spring_elementstiffness(k, theta)
    x = theta * pi / 180
    C = cos(x)
    S = sin(x)
    return k * [
        C * C C * S -C * C -C * S
        C * S S * S -C * S -S * S
        -C * C -C * S C * C C * S
        -C * S -S * S C * S S * S
    ]
end
export d2_spring_elementstiffness
"""
    function declaration: D2_TrussAssemble(K,k,i,j)

This function assembles the element stiffness
matrix k of the plane truss element with nodes
i & j into the global stiffness matrix K.
This function returns the global stiffness
matrix K after the element stiffness matrix
k is assembled.
"""
function d2_truss_assemble(K, k, i, j)
    K[2*i-1, 2*i-1] = K[2*i-1, 2*i-1] + k[1, 1]
    K[2*i-1, 2*i] = K[2*i-1, 2*i] + k[1, 2]
    K[2*i-1, 2*j-1] = K[2*i-1, 2*j-1] + k[1, 3]
    K[2*i-1, 2*j] = K[2*i-1, 2*j] + k[1, 4]
    K[2*i, 2*i-1] = K[2*i, 2*i-1] + k[2, 1]
    K[2*i, 2*i] = K[2*i, 2*i] + k[2, 2]
    K[2*i, 2*j-1] = K[2*i, 2*j-1] + k[2, 3]
    K[2*i, 2*j] = K[2*i, 2*j] + k[2, 4]
    K[2*j-1, 2*i-1] = K[2*j-1, 2*i-1] + k[3, 1]
    K[2*j-1, 2*i] = K[2*j-1, 2*i] + k[3, 2]
    K[2*j-1, 2*j-1] = K[2*j-1, 2*j-1] + k[3, 3]
    K[2*j-1, 2*j] = K[2*j-1, 2*j] + k[3, 4]
    K[2*j, 2*i-1] = K[2*j, 2*i-1] + k[4, 1]
    K[2*j, 2*i] = K[2*j, 2*i] + k[4, 2]
    K[2*j, 2*j-1] = K[2*j, 2*j-1] + k[4, 3]
    K[2*j, 2*j] = K[2*j, 2*j] + k[4, 4]
    return K
end
export d2_truss_assemble
"""
    function declaration: d2_truss_elementforce(E,A,L,theta,u)

This function returns the element force
given the modulus of elasticity E; the
cross-sectional area A; the length L;
the angle theta [in degrees], & the
element nodal displacement vector u.
"""
function d2_truss_elementforce(E, A, L, theta, u)
    x = theta * pi / 180
    C = cos(x)
    S = sin(x)
    return E * A / L * [-C -S C S] * u
end
export d2_truss_elementforce
"""
    function declaration: d2_truss_elementlength(x1,y1,x2,y2)

This function returns the length of the
plane truss element whose first node has
coordinates [x1,y1] & second node has
coordinates [x2,y2].
"""
function d2_truss_elementlength(x1, y1, x2, y2)
    return sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1))
end
export d2_truss_elementlength
"""
    function declaration: d2_truss_elementstiffness(E,A,L,theta)

This function returns the element
stiffness matrix for a plane truss
element with modulus of elasticity E;
cross-sectional area A; length L; &
angle theta [in degrees].
The size of the element stiffness
matrix is 4 x 4.
"""
function d2_truss_elementstiffness(E, A, L, theta)
    x = theta * pi / 180
    C = cos(x)
    S = sin(x)
    return E * A / L * [
        C * C C * S -C * C -C * S
        C * S S * S -C * S -S * S
        -C * C -C * S C * C C * S
        -C * S -S * S C * S S * S
    ]
end
export d2_truss_elementstiffness
"""
    function declaration: d2_truss_elementstrain(L,theta,u)

This function returns the element strain
given the length L,the angle theta [in degrees]
and the element nodal displacement vector u.
"""
function d2_truss_elementstrain(L, theta, u)
    x = theta * pi / 180
    C = cos(x)
    S = sin(x)
    return 1 / L * [-C -S C S] * u
end
export d2_truss_elementstrain
"""
    function declaration: d2_truss_elementstress(E,L,theta,u)

This function returns the element stress
given the modulus of elasticity E; the
the length L; the angle theta (in
degrees); & the element nodal
displacement vector u.
"""
function d2_truss_elementstress(E, L, theta, u)
    x = theta * pi / 180
    C = cos(x)
    S = sin(x)
    return E / L * [-C -S C S] * u
end
export d2_truss_elementstress
"""
    function declaration: d3_spring_assemble(K,k,i,j)

This function assembles the element stiffness
matrix k of the 3D spring element with nodes
i & j into the global stiffness matrix K.
This function returns the global stiffness
matrix K after the element stiffness matrix
k is assembled.
"""
function d3_spring_assemble(K, k, i, j)
    K[3*i-2, 3*i-2] = K[3*i-2, 3*i-2] + k[1, 1]
    K[3*i-2, 3*i-1] = K[3*i-2, 3*i-1] + k[1, 2]
    K[3*i-2, 3*i] = K[3*i-2, 3*i] + k[1, 3]
    K[3*i-2, 3*j-2] = K[3*i-2, 3*j-2] + k[1, 4]
    K[3*i-2, 3*j-1] = K[3*i-2, 3*j-1] + k[1, 5]
    K[3*i-2, 3*j] = K[3*i-2, 3*j] + k[1, 6]
    K[3*i-1, 3*i-2] = K[3*i-1, 3*i-2] + k[2, 1]
    K[3*i-1, 3*i-1] = K[3*i-1, 3*i-1] + k[2, 2]
    K[3*i-1, 3*i] = K[3*i-1, 3*i] + k[2, 3]
    K[3*i-1, 3*j-2] = K[3*i-1, 3*j-2] + k[2, 4]
    K[3*i-1, 3*j-1] = K[3*i-1, 3*j-1] + k[2, 5]
    K[3*i-1, 3*j] = K[3*i-1, 3*j] + k[2, 6]
    K[3*i, 3*i-2] = K[3*i, 3*i-2] + k[3, 1]
    K[3*i, 3*i-1] = K[3*i, 3*i-1] + k[3, 2]
    K[3*i, 3*i] = K[3*i, 3*i] + k[3, 3]
    K[3*i, 3*j-2] = K[3*i, 3*j-2] + k[3, 4]
    K[3*i, 3*j-1] = K[3*i, 3*j-1] + k[3, 5]
    K[3*i, 3*j] = K[3*i, 3*j] + k[3, 6]
    K[3*j-2, 3*i-2] = K[3*j-2, 3*i-2] + k[4, 1]
    K[3*j-2, 3*i-1] = K[3*j-2, 3*i-1] + k[4, 2]
    K[3*j-2, 3*i] = K[3*j-2, 3*i] + k[4, 3]
    K[3*j-2, 3*j-2] = K[3*j-2, 3*j-2] + k[4, 4]
    K[3*j-2, 3*j-1] = K[3*j-2, 3*j-1] + k[4, 5]
    K[3*j-2, 3*j] = K[3*j-2, 3*j] + k[4, 6]
    K[3*j-1, 3*i-2] = K[3*j-1, 3*i-2] + k[5, 1]
    K[3*j-1, 3*i-1] = K[3*j-1, 3*i-1] + k[5, 2]
    K[3*j-1, 3*i] = K[3*j-1, 3*i] + k[5, 3]
    K[3*j-1, 3*j-2] = K[3*j-1, 3*j-2] + k[5, 4]
    K[3*j-1, 3*j-1] = K[3*j-1, 3*j-1] + k[5, 5]
    K[3*j-1, 3*j] = K[3*j-1, 3*j] + k[5, 6]
    K[3*j, 3*i-2] = K[3*j, 3*i-2] + k[6, 1]
    K[3*j, 3*i-1] = K[3*j, 3*i-1] + k[6, 2]
    K[3*j, 3*i] = K[3*j, 3*i] + k[6, 3]
    K[3*j, 3*j-2] = K[3*j, 3*j-2] + k[6, 4]
    K[3*j, 3*j-1] = K[3*j, 3*j-1] + k[6, 5]
    K[3*j, 3*j] = K[3*j, 3*j] + k[6, 6]
    return K
end
export d3_spring_assemble
"""
    function declaration: d3_spring_elementforce(k,thetax,thetay,thetaz,u)

This function returns the element force
given the stiffness k;
the angles thetax; thetay; thetaz
(in degrees), & the element nodal
displacement vector u.
"""
function d3_spring_elementforce(k, thetax, thetay, thetaz, u)
    x = thetax * pi / 180
    w = thetay * pi / 180
    v = thetaz * pi / 180
    Cx = cos(x)
    Cy = cos(w)
    Cz = cos(v)
    return k * [-Cx -Cy -Cz Cx Cy Cz] * u
end
export d3_spring_elementforce
"""
    function declaration: d3_spring_elementstiffness(k,thetax,thetay,thetaz)

This function returns the element
stiffness matrix for a 3D spring
element with stiffness k;
angles thetax; thetay; thetaz
(in degrees). The size of the element
stiffness matrix is 6 x 6.
"""
function d3_spring_elementstiffness(k, thetax, thetay, thetaz)
    x = thetax * pi / 180
    u = thetay * pi / 180
    v = thetaz * pi / 180
    Cx = cos(x)
    Cy = cos(u)
    Cz = cos(v)
    w = [
        Cx * Cx Cx * Cy Cx * Cz
        Cy * Cx Cy * Cy Cy * Cz
        Cz * Cx Cz * Cy Cz * Cz
    ]
    return k * [w -w; -w w]
end
export d3_spring_elementstiffness
"""
    function declaration: d3_truss_assemble(K,k,i,j)

This function assembles the element stiffness
matrix k of the space truss element with nodes
i & j into the global stiffness matrix K.
This function returns the global stiffness
matrix K after the element stiffness matrix
k is assembled.
"""
function d3_truss_assemble(K, k, i, j)
    K[3*i-2, 3*i-2] = K[3*i-2, 3*i-2] + k[1, 1]
    K[3*i-2, 3*i-1] = K[3*i-2, 3*i-1] + k[1, 2]
    K[3*i-2, 3*i] = K[3*i-2, 3*i] + k[1, 3]
    K[3*i-2, 3*j-2] = K[3*i-2, 3*j-2] + k[1, 4]
    K[3*i-2, 3*j-1] = K[3*i-2, 3*j-1] + k[1, 5]
    K[3*i-2, 3*j] = K[3*i-2, 3*j] + k[1, 6]
    K[3*i-1, 3*i-2] = K[3*i-1, 3*i-2] + k[2, 1]
    K[3*i-1, 3*i-1] = K[3*i-1, 3*i-1] + k[2, 2]
    K[3*i-1, 3*i] = K[3*i-1, 3*i] + k[2, 3]
    K[3*i-1, 3*j-2] = K[3*i-1, 3*j-2] + k[2, 4]
    K[3*i-1, 3*j-1] = K[3*i-1, 3*j-1] + k[2, 5]
    K[3*i-1, 3*j] = K[3*i-1, 3*j] + k[2, 6]
    K[3*i, 3*i-2] = K[3*i, 3*i-2] + k[3, 1]
    K[3*i, 3*i-1] = K[3*i, 3*i-1] + k[3, 2]
    K[3*i, 3*i] = K[3*i, 3*i] + k[3, 3]
    K[3*i, 3*j-2] = K[3*i, 3*j-2] + k[3, 4]
    K[3*i, 3*j-1] = K[3*i, 3*j-1] + k[3, 5]
    K[3*i, 3*j] = K[3*i, 3*j] + k[3, 6]
    K[3*j-2, 3*i-2] = K[3*j-2, 3*i-2] + k[4, 1]
    K[3*j-2, 3*i-1] = K[3*j-2, 3*i-1] + k[4, 2]
    K[3*j-2, 3*i] = K[3*j-2, 3*i] + k[4, 3]
    K[3*j-2, 3*j-2] = K[3*j-2, 3*j-2] + k[4, 4]
    K[3*j-2, 3*j-1] = K[3*j-2, 3*j-1] + k[4, 5]
    K[3*j-2, 3*j] = K[3*j-2, 3*j] + k[4, 6]
    K[3*j-1, 3*i-2] = K[3*j-1, 3*i-2] + k[5, 1]
    K[3*j-1, 3*i-1] = K[3*j-1, 3*i-1] + k[5, 2]
    K[3*j-1, 3*i] = K[3*j-1, 3*i] + k[5, 3]
    K[3*j-1, 3*j-2] = K[3*j-1, 3*j-2] + k[5, 4]
    K[3*j-1, 3*j-1] = K[3*j-1, 3*j-1] + k[5, 5]
    K[3*j-1, 3*j] = K[3*j-1, 3*j] + k[5, 6]
    K[3*j, 3*i-2] = K[3*j, 3*i-2] + k[6, 1]
    K[3*j, 3*i-1] = K[3*j, 3*i-1] + k[6, 2]
    K[3*j, 3*i] = K[3*j, 3*i] + k[6, 3]
    K[3*j, 3*j-2] = K[3*j, 3*j-2] + k[6, 4]
    K[3*j, 3*j-1] = K[3*j, 3*j-1] + k[6, 5]
    K[3*j, 3*j] = K[3*j, 3*j] + k[6, 6]
    return K
end
export d3_truss_assemble
"""
    function declaration: d3_truss_elementforce(E,A,L,thetax,thetay,thetaz,u)

This function returns the element force
given the modulus of elasticity E; the
cross-sectional area A; the length L;
the angles thetax; thetay; thetaz
(in degrees), & the element nodal
displacement vector u.
"""
function d3_truss_elementforce(E, A, L, thetax, thetay, thetaz, u)
    x = thetax * pi / 180
    w = thetay * pi / 180
    v = thetaz * pi / 180
    Cx = cos(x)
    Cy = cos(w)
    Cz = cos(v)
    return E * A / L * [-Cx -Cy -Cz Cx Cy Cz] * u
end
export d3_truss_elementforce
"""
    function declaration: d3_truss_elementlength(x1,y1,z1,x2,y2,z2)

This function returns the length of the
space truss element whose first node has
coordinates [x1,y1,z1] & second node has
coordinates [x2,y2,z2].
"""
function d3_truss_elementlength(x1, y1, z1, x2, y2, z2)
    return sqrt(
        (x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1) + (z2 - z1) * (z2 - z1),
    )
end
export d3_truss_elementlength
"""
    function declaration: d3_truss_elementstiffness(E,A,L,thetax,thetay,thetaz)

This function returns the element
stiffness matrix for a space truss
element with modulus of elasticity E;
cross-sectional area A; length L; &
angles thetax; thetay; thetaz
(in degrees). The size of the element
stiffness matrix is 6 x 6.
"""
function d3_truss_elementstiffness(E, A, L, thetax, thetay, thetaz)
    x = thetax * pi / 180
    u = thetay * pi / 180
    v = thetaz * pi / 180
    Cx = cos(x)
    Cy = cos(u)
    Cz = cos(v)
    w = [
        Cx * Cx Cx * Cy Cx * Cz
        Cy * Cx Cy * Cy Cy * Cz
        Cz * Cx Cz * Cy Cz * Cz
    ]
    return E * A / L * [w -w; -w w]
end
export d3_truss_elementstiffness
"""
        function declaration: export d3_truss_elementstrain(L,thetax,thetay,thetaz,u)

This function returns the element strain
the length L; the angles thetax; thetay;
thetaz [in degrees], & the element
nodal displacement vector u.
"""
function d3_truss_elementstrain(L, thetax, thetay, thetaz, u)
    x = thetax * pi / 180
    w = thetay * pi / 180
    v = thetaz * pi / 180
    Cx = cos(x)
    Cy = cos(w)
    Cz = cos(v)
    return 1 / L * [-Cx -Cy -Cz Cx Cy Cz] * u
end
export d3_truss_elementstrain
"""
    function declaration: d3_truss_elementstress(E,L,thetax,thetay,thetaz,u)

This function returns the element stress
given the modulus of elasticity E; the
length L; the angles thetax; thetay;
thetaz [in degrees], & the element
nodal displacement vector u.
"""
function d3_truss_elementstress(E, L, thetax, thetay, thetaz, u)
    x = thetax * pi / 180
    w = thetay * pi / 180
    v = thetaz * pi / 180
    Cx = cos(x)
    Cy = cos(w)
    Cz = cos(v)
    return E / L * [-Cx -Cy -Cz Cx Cy Cz] * u
end
export d3_truss_elementstress
end #module


text:
ENV["PYTHON"] = "python3"
import PyCall
pa = PyCall.pyimport("pyarrow")
include(joinpath(dirname(pathof(Arrow)), "../test/testtables.jl"))

for (nm, t, writekw, readkw, extratests) in testtables
    nm == "unions" && continue
    println("pyarrow roundtrip: $nm")
    io = IOBuffer()
    Arrow.write(io, t; writekw...)
    seekstart(io)
    buf = PyCall.pybytes(take!(io))
    reader = pa.ipc.open_stream(buf)
    sink = pa.BufferOutputStream()
    writer = pa.ipc.new_stream(sink, reader.schema)
    for batch in reader
        writer.write_batch(batch)
    end
    writer.close()
    buf = sink.getvalue()
    jbytes = copy(reinterpret(UInt8, buf))
    tt = Arrow.Table(jbytes)
    println("pyarrow roundtrip w/ compression: $nm")
    io = IOBuffer()
    Arrow.write(io, t; compress=((:lz4, :zstd)[rand(1:2)]), writekw...)
    seekstart(io)
    buf = PyCall.pybytes(take!(io))
    reader = pa.ipc.open_stream(buf)
    sink = pa.BufferOutputStream()
    writer = pa.ipc.new_stream(sink, reader.schema)
    for batch in reader
        writer.write_batch(batch)
    end
    writer.close()
    buf = sink.getvalue()
    jbytes = copy(reinterpret(UInt8, buf))
    tt = Arrow.Table(jbytes)
end

f1 = pa.field("f1", pa.float64(), true)
f2 = pa.field("f2", pa.int64(), false)
fu = pa.field("col1", pa.union([f1, f2], "dense"))
sch = pa.schema([fu])

xs = pa.array([2.0, 4.0, PyCall.pynothing[]], type=pa.float64())
ys = pa.array([1, 3], type=pa.int64())
types = pa.array([0, 1, 0, 1, 1], type=pa.int8())
offsets = pa.array([0, 0, 1, 1, 2], type=pa.int32())
union_arr = pa.UnionArray.from_dense(types, offsets, [xs, ys])
data = [union_arr]
batch = pa.record_batch(data, names=["col1"])
sink = pa.BufferOutputStream()
writer = pa.ipc.new_stream(sink, batch.schema)
writer.write_batch(batch)
writer.close()
buf = sink.getvalue()
jbytes = copy(reinterpret(UInt8, buf))
tt = Arrow.Table(jbytes)

max_stars_repo_name: StefanKarpinski/arrow

max_issues_repo_name: StefanKarpinski/arrow

max_forks_repo_name: StefanKarpinski/arrow


text:
using SafeTestsets

@safetestset "My f tests" begin include("my_f_tests.jl") end


text:
const to = TimerOutput()

function run( nstep; npm = 100 )

    dt = 0.01
    alpha = 0.1
    kx   = 0.5
    dimx = 2*pi/kx
    dimy = 1
    nx   = 128  # nombre de pts suivant x
    ny   = 16   # nombre de pts suivant y
    mesh = TwoDGrid( dimx, nx, dimy, ny)
    dx, dy = mesh.dx, mesh.dy

    nbpart = npm*nx*ny
    println( " nbpart = $nbpart ")

    particles = ParticleGroup{2,2}(nbpart, charge = 1.0, mass = 1.0, n_weights = 1)

    sampler = LandauDamping(alpha, kx)

    sample!(particles, mesh, sampler)

    fdtd = FDTD(mesh)
    for i=1:nx, j=1:ny+1
        fdtd.ex[i,j] = alpha/kx * sin(kx*(mesh.x[i]+mesh.x[i+1])/2)
    end
    time = 0
    energy = Float64[compute_energy(fdtd, mesh)]
    t = Float64[time]

    kernel = CloudInCell()

    reset_timer!()
    @showprogress 1 for istep in 1:nstep

       if istep > 1
           @timeit to "fdtd" faraday!( fdtd, mesh, 0.5dt )
       end
       update_fields!(mesh, fdtd)
       #@timeit to "interpolation" interpolation!(particles, mesh)
       @timeit to "pushv" push_v!( particles, kernel, mesh, dt )
       @timeit to "pushx" push_x!( particles, mesh, 0.5dt)
       @timeit to "deposition" compute_current!( mesh, kernel, particles)
       @timeit to "pushx" push_x!( particles, mesh, 0.5dt)
       @timeit to "fdtd" faraday!(fdtd, mesh, 0.5dt)
       @timeit to "fdtd" ampere_maxwell!(fdtd, mesh, dt)
       time = time + dt
       push!(t, time)
       push!(energy, compute_energy(fdtd, mesh))

    end

    t, energy

end


text:
using SGDBabysitter
using Test

@testset "SGDBabysitter.jl" begin
    # Write your own tests here.
end


text:
# SPDX-License-Identifier: MIT

# Constant parameter.
#
# Lucas Ondel 2021

"""
    struct ConstParam{T<:AbstractVector} <: AbstractParameter
        μ::T
    end

Constant parameter.

# Constructor

    ConstantParam(μ)

Create a constant parameter.
"""
struct ConstParameter{T<:AbstractVector} <: AbstractParameter
    μ::T
    _stats::Function
end

statistics(param::ConstParameter) = param._stats(param.μ)

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
"""
Find optimum reparam between two images

    reparam_image(It::Array,Im::Array,gam::Array,b::Array;
                  stepsize::Float64=1e-5, itermax::Int32=1000,
                  lmark::Bool=false)
    :param It: Template Image
    :param Im: Test Image
    :param gam: initial warping
    :param b: basis matrix
    :param stepsize: gradient stepsize
    :param itermax: maximum number of iterations
    :param lmark: use landmarks

    :return gamnew: final warping
    :return Inew: parameterized image
    :return H: energy
    :return stepsize: final stepsize
"""
function reparam_image(It::Array, Im::Array, gam::Array, b::Array;
                       stepsize::Float64=1e-5, itermax::Integer=1000,
                       lmark::Bool=false)
    m = size(It,1);
    n = size(It,2);
    gamid = makediffeoid(m,n);

    # Main Loop
    H = zeros(itermax+1);
    Iold = apply_gam_to_imag(Im,gam);
    Iold -= minimum(Iold);
    Iold /= maximum(Iold);
    gamold = copy(gam);
    qt = imag_to_q(It);
    qm = imag_to_q(Iold);

    gamnew = copy(gamold);
    Inew = copy(Iold);
    iter = 1;
    H[iter] = comp_energy(qt,qm);
    @printf("Iteration %d, energy %f\n",iter-1,H[iter])

    gamupdate = update_gam(qt,qm,b);
    cutoff = 1e-3;
    for iter = 2:(itermax+1)
        gaminc = gamid + stepsize.*gamupdate;
        G = check_crossing(gamnew);
        if !G
            @printf("Possible Crossing!")
            gamnew = copy(gamold);
            stepsize = 0.67*stepsize;
            H[iter] = H[iter-1];
            continue
        else
            gamnew = apply_gam_to_gam(gamnew, gaminc);
        end

        Inew = apply_gam_to_imag(Im, gamnew);
        Inew -= minimum(Inew);
        Inew /= maximum(Inew);
        qm = imag_to_q(Inew);
        H[iter] = comp_energy(qt,qm);
        @printf("Iteration %d, energy %f\n",iter-1,H[iter])

        if (iter > 4)
            hstop = 1.;
            for i = 1:4
                hstop = hstop*(H[iter] >= H[iter-1]);
            end
            if (hstop!=0)
                @printf("Warning: energy constantly increasing")
                break
            end
        end

        if (iter > 4) && (H[iter]>=H[iter-1]) && (H[iter-1]>=H[iter-2]) &&
            (H[iter-2] >= H[iter-3])
            @printf("Warning: energy is not changing")
        end

        if ((iter>1) && (H[iter]>H[iter-1])) || ((iter>3) &&
            ((H[iter-1]<=H[iter-2]) && (H[iter-2]>H[iter-3])))
            stepsize = 0.9*stepsize;
            gamnew = copy(gamold);
            H[iter] = H[iter-1];
            continue
        end

        gamold = copy(gamnew);
        gamupdate = update_gam(qt, qm, b);
    end

    H = H[1:iter];

    return gamnew, Inew, H, stepsize
end

"""
Pairwise align two images

    pair_align_image(I1, I2; M=5, ortho=true, basis_type="t", resizei=true,
                     N=64, stepsize=1e-5, itermax=1000)
    :param I1: reference image
    :param I2: image to warp
    :param M: number of basis elements
    :param ortho: orthonormalize basis
    :param basis_type: type of basis ("t", "s", "i", "o")
    :param resizei: resize image
    :param N: size of resized image
    :param stepsize: gradient stepsize
    :param itermax: maximum number of iterations

    :return I2_new: aligned I2
    :return gam: warping function
"""
function pair_align_image(I1, I2; M=5, ortho=true, basis_type="t", resizei=true,
                          N=64, stepsize=1e-5, itermax=1000)
    m,n = size(I1);
    F1 = zeros(m,n,2);
    m1,n1 = size(I2);
    F2 = zeros(m,n,2);

    # Take Gradient
    fu, fv = gradient2(I1,1./(m-1),1./(n-1));
    F1[:,:,1] = fu;
    F1[:,:,2] = fv;
    fu, fv = gradient2(I2,1./(m1-1),1./(n1-1));
    F2[:,:,1] = fu;
    F2[:,:,2] = fv;

    # Resize Data and Center
    if resizei
        if N > m || N > n
            @printf("Not resizing, N is larger than image size")
        else
            m_n = LinRange(1,m,N);
            n_n = LinRange(1,n,N);
            F1a = zeros(N,N,2);
            spl = Spline2D(collect(1:m),collect(1:n),F1[:,:,1]);
            F1a[:,:,1] = evalgrid(spl,m_n,n_n);
            spl = Spline2D(collect(1:m),collect(1:n),F1[:,:,2]);
            F1a[:,:,2] = evalgrid(spl,m_n,n_n);
            F1 = copy(F1a);

            m_n = LinRange(1,m1,N);
            n_n = LinRange(1,n1,N);
            F2a = zeros(N,N,2);
            spl = Spline2D(collect(1:m1),collect(1:n1),F2[:,:,1]);
            F2a[:,:,1] = evalgrid(spl,m_n,n_n);
            spl = Spline2D(collect(1:m1),collect(1:n1),F2[:,:,2]);
            F2a[:,:,2] = evalgrid(spl,m_n,n_n);
            F2 = copy(F2a);
        end
    end
    F1 -= minimum(F1);
    F1 /= maximum(F1);
    F2 -= minimum(F2);
    F2 /= maximum(F2);

    # Generate basis
    b, gamid = run_basis(F1, M, basis_type, ortho);
    gamp = copy(gamid);
    gam, F2_new, H, stepsize = reparam_image(F1, F2, gamp, b,
                                             stepsize=stepsize,
                                             itermax=itermax);

    I2_new = apply_gam_to_imag(I2, gam);

    return I2_new, gam
end


text:
function add_data_points!(df_target, ds::Dataset, p)
    for tsi in ds.timeseries
        add_data_points!(df_target,tsi, p)
    end
    for pdi in ds.ponctual
        add_data_points!(df_target,pdi, p)
    end
    return df_target
end

function add_data_points!(df_target, ts::TimeseriesDataset, p)
    if hasproperty(p.estimate,:extra_misfit_target)
        for extra_misfit_target_name in propertynames(p.estimate.extra_misfit_target)
            if extra_misfit_target_name ∈ ts.target
                dfg = group(ts)
                for g in dfg
                    extra_data_props = getproperty(p.estimate.extra_misfit_target,extra_misfit_target_name)
                    for i in eachindex(extra_data_props)
                        func = extra_data_props[i].f
                        colname = extra_misfit_target_name
                        val_target = func(g[!,colname])
                        val_target_sim = func(g[!,string(colname)*"_sim"])
                        for iv in eachindex(val_target)
                            row = (
                                target = val_target[iv],
                                target_sim = val_target_sim[iv],
                                std = extra_data_props[i].std,
                                target_name = extra_data_props[i].transform_target_name
                            )
                            push!(df_target,row)
                        end
                    end
                end
            end
        end
    end
    return df_target
end

function add_data_points!(df_target, pd::PonctualDataset, p)
    data = pd.data
    if hasproperty(p.estimate,:extra_misfit_target)
        for extra_misfit_target_name in propertynames(p.estimate.extra_misfit_target)
            if extra_misfit_target_name ∈ pd.target
                extra_data_props = getproperty(p.estimate.extra_misfit_target,extra_misfit_target_name)
                for i in eachindex(extra_data_props)
                    func = extra_data_props[i].f
                    colname = extra_misfit_target_name
                    val_target = func(data[!,colname])
                    val_target_sim = func(data[!,string(colname)*"_sim"])
                    for iv in eachindex(val_target)
                        row = (
                            target = val_target[iv],
                            target_sim = val_target_sim[iv],
                            std = extra_data_props[i].std,
                            target_name = extra_data_props[i].transform_target_name
                        )
                        push!(df_target,row)
                    end
                end
            end
        end
    end
    return df_target
end

max_stars_repo_name: Leooop/ParametersEstimator

max_issues_repo_name: Leooop/ParametersEstimator

max_forks_repo_name: Leooop/ParametersEstimator


text:
# tests for distributions defined in VinDsl
using FactCheck
using Distributions
using PDMats

using VinDsl
srand(12349)

facts("Checking MatrixNormal distribution") do
    n, p = 5, 7
    LU = rand(n, n)
    LV = rand(p, p)
    U = PDMat(LU * LU')
    V = PDMat(LV * LV')
    M = rand(n, p)

    X = MatrixNormal(M, U, V)

    context("Constructors: 3-param") do
        @fact isa(X, MatrixNormal) --> true
        @fact isa(X, ContinuousMatrixDistribution) --> true
        @fact isa(X, Distribution) --> true
        @fact X.M --> M
        @fact X.U --> U
        @fact X.V --> V
    end

    context("Constructors: 2-param, same type") do
        XX = MatrixNormal(U, V)
        @fact isa(XX, MatrixNormal) --> true
        @fact XX.M --> zeros(dim(U), dim(V))
    end

    context("Constructors: 2-param, different type") do
        UU = ScalMat(n, 2.4)
        XX = MatrixNormal(U, V)
        @fact isa(XX, MatrixNormal) --> true
        @fact XX.M --> zeros(dim(UU), dim(V))
    end

    context("Constructors: 2-param, matrix and vector") do
        XX = MatrixNormal(U.mat, V.mat)
        @fact isa(XX, MatrixNormal) --> true

        XX = MatrixNormal(U, V.mat)
        @fact isa(XX, MatrixNormal) --> true

        XX = MatrixNormal(U.mat, V)
        @fact isa(XX, MatrixNormal) --> true

        XX = MatrixNormal(diag(U.mat), V)
        @fact isa(XX, MatrixNormal) --> true

        XX = MatrixNormal(diag(U.mat), V.mat)
        @fact isa(XX, MatrixNormal) --> true
    end

    context("Bad constructor inputs throw errors") do
        @fact_throws DimensionMismatch MatrixNormal(M[2:end, :], U, V)
        @fact_throws DimensionMismatch MatrixNormal(M[:, 2:end], U, V)
    end

    context("Check basic interface") do
        @fact size(X) --> (n, p)
        @fact params(X) --> (M, U, V)
    end

    context("Check mean") do
        m = mean(X)
        @fact m --> X.M
    end

    context("Check cov and var") do
        @fact cov(X) --> kron(V.mat, U.mat)
        @fact var(X) --> diag(U.mat) * diag(V.mat)'
        @fact invcov(X) --> roughly(invcov(convert(MvNormal, X)))
        @fact logdetcov(X) --> roughly(logdet(kron(full(V), full(U))))
    end

    context("Check sampling") do
        z = rand(X)
        @fact size(z) --> size(X)
    end

    context("Check logpdf") do
        z = rand(X)
        lpdf = logpdf(X, z)
        m = convert(MvNormal, X)

        @fact lpdf --> roughly(logpdf(m, vec(z)))
    end

    context("Check entropy") do
        @fact entropy(X) --> entropy(convert(MvNormal, X))
    end

end

facts("Checking exponential family interface") do
    context("Normal") do
        μ = 1.1
        σ = 2.5
        τ = 1/σ^2
        d = Normal(μ, σ)
        @fact naturals(d) --> (μ * τ, -τ/2)
        @fact naturals_to_params(naturals(d), Normal) --> (μ, σ)
        @fact Normal(constrain(unconstrain(d), d)...) --> d
        @fact supp(d) --> RReal()
    end

    context("Gamma") do
        a = 1.1
        θ = 2.5
        d = Gamma(a, θ)
        @fact naturals(d) --> (a - 1, -1/θ)
        @fact naturals_to_params(naturals(d), Gamma) --> (a, θ)
        @fact Gamma(constrain(unconstrain(d), d)...) --> d
        @fact supp(d) --> RPositive()
    end

    context("Dirichlet") do
        a = rand(5)
        α = a / sum(a)
        d = Dirichlet(α)
        @fact naturals(d) --> (α - 1,)
        @fact naturals_to_params(naturals(d), Dirichlet)[1] --> roughly(α)
        # @fact Dirichlet(constrain(unconstrain(d), Dirichlet)...) --> d
    end

    context("MvNormalCanon") do
        p = 7
        h = rand(p)
        jj = rand(p, p)
        J = jj * jj'
        d = MvNormalCanon(h, J)
        @fact naturals(d) --> (h, -J/2)
        @fact naturals_to_params(naturals(d), MvNormalCanon) --> (h, J)
        dd =  MvNormalCanon(constrain(unconstrain(d), d)...)
        @fact dd.h --> h
        @fact dd.J.mat --> roughly(J)
        @fact supp(d) --> RRealVec(length(d))
    end

    context("Wishart") do
        p = 7
        ss = rand(p, p)
        S = ss * ss'
        df = p + 10
        d = Wishart(df, S)
        nats = naturals(d)
        @fact nats[1] --> (df - p - 1)/2
        @fact nats[2] --> roughly(-inv(S)/2)
        n2p =  naturals_to_params(naturals(d), Wishart)
        @fact n2p[1] --> df
        @fact n2p[2] --> roughly(S)
        dd = Wishart(constrain(unconstrain(d), d)...)
        @fact dd.df --> df
        @fact dd.S.mat --> roughly(S)
        @fact supp(d) --> RCovMat(dim(d))
    end
end

facts("Checking LKJ distribution") do
    lkjeta = 4
    p = 5
    yy = LKJ(lkjeta, p)
    @fact dim(yy) --> 5
    @fact size(yy) --> (5, 5)
    @fact params(yy) --> (4, 5)
    @fact typeof(yy) --> LKJ{Float64}
    @fact partype(yy) --> Float64
    @fact meanlogdet(yy) --> -7.462607436560608 # double checked with John's Python code after minor fixes
    @fact meanlog(yy, 10) --> 2.5154263517968545
    @fact entropy(yy) --> -4.787291422054154
end


text:
export overlay_mpact_grid_hierarchy

function overlay_mpact_grid_hierarchy(grid::MPACTGridHierarchy,
                                      material_hierarchy::Vector{Material})
    @info "Overlaying MPACT grid hierarchy"
    # Generate all of the rectangles the make up the coarse grid, and save
    # their (x,y) origin coordinates, so we can group them into their appropriate
    # modules and lattices later
    model_dtags = gmsh.model.get_entities(2)
    lg_x = grid.lattice_grid.x
    lg_y = grid.lattice_grid.y
    lg_nx = length(lg_x)
    lg_ny = length(lg_y)
    mg_x = grid.module_grid.x
    mg_y = grid.module_grid.y
    mg_nx = length(mg_x)
    mg_ny = length(mg_y)
    cg_x = grid.coarse_grid.x
    cg_y = grid.coarse_grid.y
    cg_nx = length(cg_x)
    cg_ny = length(cg_y)
    grid_tags_coords = Vector{Tuple{Int32,Float64,Float64}}(undef, (cg_nx-1)*(cg_ny-1))
    gtctr = 1
    for iy in 1:cg_ny - 1
        y = cg_y[iy]
        for ix in 1:cg_nx - 1
            x = cg_x[ix]
            tag = gmsh.model.occ.add_rectangle(x, y, 0, cg_x[ix+1] - x, cg_y[iy+1] - y)
            grid_tags_coords[gtctr] = (tag, x, y)
            gtctr += 1
        end
    end
    gmsh.model.occ.synchronize()
    grid_dtags = [(Int32(2), gtc[1]) for gtc in grid_tags_coords]

    # Label the rectangles with the appropriate physical groups
    # Create a dictionary holding all the physical group names and tags
    groups = Dict{String, Vector{Int32}}()
    max_grid_digits = max(length(string(cg_nx-1)),
                          length(string(cg_ny-1)))

    # Create each grid name
    # Lattices
    for ix in 1:lg_nx - 1
        for iy in 1:lg_ny - 1
            grid_str = string("Lattice (", lpad(ix, max_grid_digits, "0"), ", ",
                                           lpad(iy, max_grid_digits, "0"), ")")
            groups[grid_str] = Int32[]
        end
    end

    # Modules
    for ix in 1:mg_nx - 1
        for iy in 1:mg_ny - 1
            grid_str = string("Module (", lpad(ix, max_grid_digits, "0"), ", ",
                                          lpad(iy, max_grid_digits, "0"), ")")
            groups[grid_str] = Int32[]
        end
    end

    # Coarse Cells
    for ix in 1:cg_nx - 1
        for iy in 1:cg_ny - 1
            grid_str = string("Coarse Cell (", lpad(ix, max_grid_digits, "0"), ", ",
                                               lpad(iy, max_grid_digits, "0"), ")")
            groups[grid_str] = Int32[]
        end
    end

    # For each rectangle, find which lattice/module/coarse cell it belongs to
    # Lattices
    for (tag, x, y) in grid_tags_coords
        ix = searchsortedfirst(lg_x, x)
        if x == lg_x[ix]
            ix += 1
        end
        iy = searchsortedfirst(lg_y, y)
        if y == lg_y[iy]
            iy += 1
        end
        grid_str = string("Lattice (", lpad(ix-1, max_grid_digits, "0"), ", ",
                                       lpad(iy-1, max_grid_digits, "0"), ")")
        push!(groups[grid_str], tag)
    end

    # Modules
    for (tag, x, y) in grid_tags_coords
        ix = searchsortedfirst(mg_x, x)
        if x == mg_x[ix]
            ix += 1
        end
        iy = searchsortedfirst(mg_y, y)
        if y == mg_y[iy]
            iy += 1
        end
        grid_str = string("Module (", lpad(ix-1, max_grid_digits, "0"), ", ",
                                      lpad(iy-1, max_grid_digits, "0"), ")")
        push!(groups[grid_str], tag)
    end

    # Coarse cells
    for (tag, x, y) in grid_tags_coords
        ix = searchsortedfirst(cg_x, x)
        if x == cg_x[ix]
            ix += 1
        end
        iy = searchsortedfirst(cg_y, y)
        if y == cg_y[iy]
            iy += 1
        end
        grid_str = string("Coarse Cell (", lpad(ix-1, max_grid_digits, "0"), ", ",
                                           lpad(iy-1, max_grid_digits, "0"), ")")
        push!(groups[grid_str], tag)
    end

    # Setup material physical group
    mat_name = "Material: "*material_hierarchy[end].name
    groups[mat_name] = [dt[2] for dt in grid_dtags]
    old_groups = gmsh.model.get_physical_groups()
    names = [gmsh.model.get_physical_name(grp[1], grp[2]) for grp in old_groups]

    # Assign groups
    for name in keys(groups)
        safe_add_physical_group(name, [(Int32(2),t) for t in groups[name]])
    end

    # material hierarchy with the grid material at the bottom.
    output_dtags, output_dtags_map = safe_fragment(
                        model_dtags, grid_dtags,
                        material_hierarchy = material_hierarchy
                   )

    # Account for material already in model, and erase an add.
    return output_dtags, output_dtags_map
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using KernelFunctions: KernelFunctions
using PyCall: PyCall

# use Julia wrapper for ColVecs and RowVecs in Python
# otherwise the automatic conversions in PyCall will convert it to a matrix!
function PyCall.PyObject(::Type{<:KernelFunctions.ColVecs})
    return PyCall.pyfunctionret(KernelFunctions.ColVecs, Any, PyCall.PyAny)
end
function PyCall.PyObject(::Type{<:KernelFunctions.RowVecs})
    return PyCall.pyfunctionret(KernelFunctions.RowVecs, Any, PyCall.PyAny)
end

max_stars_repo_name: devmotion/pycalibration

max_issues_repo_name: devmotion/pycalibration

max_forks_repo_name: devmotion/pycalibration