text:
function is_valid(reactor, R_0, B_0)
  isnan(R_0) && return false

  ( R_0 < R_0_min ) && return false
  ( R_0 > R_0_max ) && return false

  checked_fields = [:H, :I_M, :n_20, :nu_J, :T_k, :B_C]

  tmp_is_valid = true

  for checked_field in checked_fields
    work_value = getfield(
      FusionSystems, Symbol("calc_$(checked_field)")
    )(reactor, R_0, B_0)

    work_max = getfield(
      FusionSystems, Symbol("$(checked_field)_max")
    )

    if isdefined(FusionSystems, Symbol("$(checked_field)_min"))
      work_min = getfield(FusionSystems, Symbol("$(checked_field)_min"))
    else
      work_min = checked_min
    end

    ( work_min < work_value < work_max ) && continue
    tmp_is_valid = false
    break
  end

  tmp_is_valid
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualNetworkPropertiesFormat <: SwaggerModel
    addressSpace::Any # spec type: Union{ Nothing, AddressSpace } # spec name: addressSpace
    dhcpOptions::Any # spec type: Union{ Nothing, DhcpOptions } # spec name: dhcpOptions
    subnets::Any # spec type: Union{ Nothing, Vector{Subnet} } # spec name: subnets
    virtualNetworkPeerings::Any # spec type: Union{ Nothing, Vector{VirtualNetworkPeering} } # spec name: virtualNetworkPeerings
    resourceGuid::Any # spec type: Union{ Nothing, String } # spec name: resourceGuid
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function VirtualNetworkPropertiesFormat(;addressSpace=nothing, dhcpOptions=nothing, subnets=nothing, virtualNetworkPeerings=nothing, resourceGuid=nothing, provisioningState=nothing)
        o = new()
        validate_property(VirtualNetworkPropertiesFormat, Symbol("addressSpace"), addressSpace)
        setfield!(o, Symbol("addressSpace"), addressSpace)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("dhcpOptions"), dhcpOptions)
        setfield!(o, Symbol("dhcpOptions"), dhcpOptions)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("subnets"), subnets)
        setfield!(o, Symbol("subnets"), subnets)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("virtualNetworkPeerings"), virtualNetworkPeerings)
        setfield!(o, Symbol("virtualNetworkPeerings"), virtualNetworkPeerings)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("resourceGuid"), resourceGuid)
        setfield!(o, Symbol("resourceGuid"), resourceGuid)
        validate_property(VirtualNetworkPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type VirtualNetworkPropertiesFormat

const _property_map_VirtualNetworkPropertiesFormat = Dict{Symbol,Symbol}(Symbol("addressSpace")=>Symbol("addressSpace"), Symbol("dhcpOptions")=>Symbol("dhcpOptions"), Symbol("subnets")=>Symbol("subnets"), Symbol("virtualNetworkPeerings")=>Symbol("virtualNetworkPeerings"), Symbol("resourceGuid")=>Symbol("resourceGuid"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_VirtualNetworkPropertiesFormat = Dict{Symbol,String}(Symbol("addressSpace")=>"AddressSpace", Symbol("dhcpOptions")=>"DhcpOptions", Symbol("subnets")=>"Vector{Subnet}", Symbol("virtualNetworkPeerings")=>"Vector{VirtualNetworkPeering}", Symbol("resourceGuid")=>"String", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ VirtualNetworkPropertiesFormat }) = collect(keys(_property_map_VirtualNetworkPropertiesFormat))
Swagger.property_type(::Type{ VirtualNetworkPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualNetworkPropertiesFormat[name]))}
Swagger.field_name(::Type{ VirtualNetworkPropertiesFormat }, property_name::Symbol) =  _property_map_VirtualNetworkPropertiesFormat[property_name]

function check_required(o::VirtualNetworkPropertiesFormat)
    true
end

function validate_property(::Type{ VirtualNetworkPropertiesFormat }, name::Symbol, val)
end

ext: jl


text:
let

install = Base.find_in_path("Atom") == nothing

if install
  print(STDERR, "juno-msg-installing")
  try
    Pkg.add("Atom")
  catch
    print(STDERR, "juno-msg-install")
    rethrow()
  end
end

port = parse(Int, shift!(ARGS))

precompile = !install && isempty(Base.find_all_in_cache_path(:Atom))

if precompile
  print(STDERR, "juno-msg-precompiling")
end

junorc = abspath(homedir(), ".junorc.jl")

try
  import Atom
  using Juno
  Atom.connect(port, welcome = precompile || install)
  Atom.handle("connected") do
    ispath(junorc) && include(junorc)
    nothing
  end
catch
  print(STDERR, "juno-msg-load")
  rethrow()
end

end

max_stars_repo_name: jessebett/atom-julia-client

max_issues_repo_name: jessebett/atom-julia-client

max_forks_repo_name: jessebett/atom-julia-client


text:
module Arcpy

using PyCall
@pyimport arcpy

include("toolbox.jl")

GetParameterAsText(x) = arcpy.GetParameterAsText(x)
GetRasterProperties_management(x, y) = arcpy.GetRasterProperties_management(x, y)
getOutput(x, i) = x[:getOutput](i)

AddMessage(str) = arcpy.AddMessage(str)
setenv(var, val) = (arcpy.env[Symbol(var)] = val)

for f in (:FeatureClassToFeatureClass_conversion,
		:Dissolve_management, :Delete_management,
		:Clip_management, :CreateFishnet_management,
		:AddField_management, :CalculateField_management,
		:AddFieldDelimiters)
	@eval ($f)(args...) = acrpy.($f)(args...)
end

end # module


text:
using Plots, SeisNoise, DTWDT, PlotlyJS
"""
Allen and Smith autoregressive lag-1 autocorrelation coefficient.
In an AR(1) model
    x(t) - <x> = γ(x(t-1) - <x>) + α z(t)
where <x> is the process mean, γ and α are process parameters and z(t) is a
Gaussian unit-variance white noise

References
    ----------
    [1] Allen, M. R. and Smith, L. A. Monte Carlo SSA: detecting
        irregular oscillations in the presence of colored noise.
        *Journal of Climate*, **1996**, 9(12), 3373-3404.
        <http://dx.doi.org/10.1175/1520-0442(1996)009<3373:MCSDIO>2.0.CO;2>
    [2] http://www.madsci.org/posts/archives/may97/864012045.Eg.r.html

Original function translated from https://github.com/regeirk/pycwt/blob/master/pycwt/helpers.py
"""
function ar1(x::AbstractArray)
    N = size(x)
    xm = mean(x)
    x = x .- xm

    # estimate the lag zero and one covariance
    c0 = (x' .* x) / N
    c1 = (x[0:N-1]' .* x[1:N]) / (N-1)

    # according to A. Grinsteds' substitutions
    A = c0 * N^2
    B = -c1 .* N .- A - 2*c0 .+ 2*c1 .- c1 * N^2 .+ c0*N
    C = N * (c0 .+ c1*N .- c1)
    D = B^2 - 4*A .* C

    if D > 0
        g = (-B - sqrt(D) / (2*A))
    else
        @warn "Cannot place an upperbound on the unbiased AR(1). Series is
               too short or trend too large."
    end

    # according to Allen & Smith (1996), footnote 4
    μ2 = -1 / N + (2 / N^2) * ((N - g^N) / (1-g) - g * (1 - g^(N-1)) / (1-g)^2)
    c0t = c0 / (1-μ2)
    a = sqrt((1-g^2) * c0t)

    return g, a, μ2
end

function smooth(W::AbstractArray, dt::Float64, dj::Float64, scales::Array{Float64,1})
    # The smoothing is performed by using a filter given by the absolute
    # value of the wavelet function at each scale, normalized to have a
    # total weight of unity, according to suggestions by Torrence &
    # Webster (1999) and by Grinsted et al. (2004).
    m, n = size(W)

    # Filter in time
    k = 2π*fftfreq(length(W[1,:]))
    k2 = k.^2
    snorm = scales ./ dt

    # Smoothing by Gaussian window (absolute value of wavelet function)
    # using the convolution theorem: multiplication by Gaussian curve in
    # Fourier domain for each scale, outer product of scale and frequency
    F = exp.(-0.5 .* (snorm.^2) .* k2') # outer product
    smooth = (ifft(F .* fft(W,2), 2))

    T = smooth[:, 1:n] # Remove possible padded region due to FFTW

    # Filter in scale. For the Morlet wavelet, this is simply a boxcar with 0.6 width
    # construct boxcar
    wsize = convert(Int64, round(0.60 / dj * 2))
    if wsize % 2 == 0
        wsize+=1
    end
    halfWin = div(wsize,2)

    # iterate over 'horizontal' and smooth in the 'vertical'
    # this could also be done by adding an axis to the transpose and performing a 2d convolution
    for i=1:size(T,2)
        # pad signal for 'same' padding
        paddedT = vcat(T[1, i]*ones(halfWin), T[:, i], T[end, i]*ones(halfWin))
        # smooth signal
        win = Array{eltype(paddedT), 1}(undef, wsize)
        win.=ones(wsize)/wsize
        smoothT = conv(paddedT, win)
        # trim signal
        T[:, i] = smoothT[2*halfWin+1:end-2*halfWin]
    end

    return T
end

"""
Modified wavelet coherence transform
Modified from original by Congcong Yuan
"""
function wct_modified(y1, y2, dt, dj, s0, J; wavelet::String="morlet", norm::Bool=true, sig::Bool=false, siglvl::Float64=0.95)
    # define wavelet
    if wavelet=="morlet"
        wav = WT.Morlet(6)
    end

    fλ = (4*π) / (wav.σ + sqrt(2 + wav.σ^2))

    # calculate std dev of both input signals (std(y1), std(y2))
    std1 = std(y1)
    std2 = std(y2)
    # normalize
    if norm
        y1 = (y1 .- mean(y1)) ./ std1
        y2 = (y2 .- mean(y2)) ./ std2
    end

    # calculate the CWT of the time series, using identical parameters for both calculations
    W1, sj, freqs, coi = cwt(y1, CFW(wav,1/dj), J1=J, dt=dt, s0=s0)
    W2, sj, freqs, coi = cwt(y2, CFW(wav,1/dj), J1=J, dt=dt, s0=s0)

    # fourier frequencies corresponding to wavelet scales
    freqs = 1 ./ (fλ .* sj)

    scales1 = Array{Float64,2}(undef, size(W1))
    scales2 = Array{Float64,2}(undef, size(W1))
    for i=1:size(scales1, 2)
        scales1[:,i] .= sj
        scales2[:,i] .= sj
    end

    # smooth wavelet spectra before truncating
    S1 = smooth(abs.(W1).^2 ./ scales1, dt, dj, sj)
    S2 = smooth(abs.(W2).^2 ./ scales2, dt, dj, sj)

    # compute wavelet transform coherence
    W12 = W1 .* conj(W2)

    scales = scales1
    S12 = smooth(W12 ./ scales, dt, dj, sj)
    WCT = abs.(S12).^2 ./ (S1 .* S2)
    aWCT = angle.(W12)

    # calculate cross spectrum and its amplitude
    WXS = W12
    WXA = abs.(S12)

    # calculate the significance using a Monte Carlo simulation with 95%
    # confidence as a function of scale
    if sig
        # TODO: implement and test the wct_significance and ar1 functions
        a1, b1, c1 = ar1(y1)
        a2, b2, c2 = ar1(y2)
        sig = wct_significance(a1, a2, dt=dt, dj=dj, s0=s0, J=J, significance_level=siglvl, wavelet=wavelet)
    else
        sig = [0]
    end

    return WXS, WXA, real.(WCT), aWCT, coi, freqs, sig
end

"""
Compute wavelet cross spectrum given an array of frequency band minima and maxima.
"""
function wxs_freqbands(cur, ref, t, twindow, freqbands, dj, s0, J; wvn::String="morlet", unwrapflag::Bool=false, sig::Bool=false, siglvl::Float64=0.95, windowflag::Bool=true, nwindow::Float64=1.5)
    # define sample frequency
    dt = t[2] - t[1]
    fs = 1/dt

    # perform wavelet coherence transform
    WXS, WXA, WCT, aWCT, coi, freqs, sig = wct_modified(cur, ref, dt, dj, s0, J, wavelet=wvn, sig=sig, siglvl=siglvl)

    # do inverse cwt for different frequency bands
    if unwrapflag==true
        phase = unwrap(aWCT, dims=ndims(aWCT))
    else
        phase = aWCT
    end

    (nband, nfreq) = size(freqbands)
    if nband==0 || nfreq!=2
        println("Error: please check if inputs in freqbands is right!")
    else
        dtt = zeros(nband)
        err = zeros(nband)

        # iterate over frequency bands
        for iband=1:nband
            fmin = freqbands[iband, 1]
            fmax = freqbands[iband, 2]
            if (fmax > maximum(freqs)) || (fmax <= fmin)
                println("Error: please ensure columns 1 and 2 are right frequency limits in freqbands!")
            else
                freq_ind = findall(f->(f>=fmin && f<=fmax), freqs)

                iphase = []
                iphase = phase[freq_ind, :]

                if windowflag
                    tmin = twindow[1]
                    tmax = twindow[2]

                    # time checks
                    if tmin < minimum(t) || tmax > maximum(t) || tmax <= tmin
                        println("Error: please input right time limits in the time window!")
                    else
                        # truncate data with the time window
                        t_ind = findall(x->(x≤tmax && x≥tmin), t)
                        wt = t[t_ind]
                    end
                else # using dynamic time window
                    tmin = twindow[1]
                    tmax = twindow[1] + nwindow*(1/fmin + 1/fmax)/2.

                    if tmin < minimum(t) || tmax > maximum(t) || tmax <= tmin
                        println("Error: please input right time limits in the time window!")
                    else
                        # truncate data with the time window
                        t_ind = findall(x->(x≤tmax && x≥tmin), t)
                        wt = t[t_ind]
                    end
                end

                delta_t = zeros(nband, length(t))
                delta_t_err = zeros(nband, length(t))
                # do regression to get slope dt_m between phase delays and frequency band
                for itime=1:length(wt)
                    # get weights
                    # w = 1 ./ WXA[freq_ind,itime]
                    w = 1 ./ WCT[freq_ind, itime]
                    infNaN = findall(x->(isnan.(x) || isinf.(x)), w)
                    if length(infNaN)!=0
                        w[infNaN] .= 1.0
                    end

                    # perform a WLS inversion
                    # This does NOT force the best fit line through the origin
                    model = glm(@formula(Y ~ X),DataFrame(X=freqs[freq_ind]*2π,Y=iphase[:,itime]),Normal(),
                    IdentityLink(),wts=w)

                    delta_t[iband, itime] = coef(model)[2]
                    delta_t_err[iband, itime] = stderror(model)[2]
                end

                # do regression in time to get time shift
                # w2 = 1/np.mean(WXA[freq_ind,:],axis=0)
                w2 = 1 ./ mean(WCT[freq_ind,:], dims=1)
                infNaN =findall(x->(isnan.(x) || isinf.(x)), w2[:])
                if length(infNaN)!=0
                    w2[:, infNaN] .= 1.0
                end
            end

            if length(wt)>2
                if delta_t == nothing
                    continue
                end

                # perform a WLS inversion
                # This DOES force the best fit line through the origin
                model = glm(@formula(Y ~0 + X),DataFrame(X=wt,Y=delta_t[iband, :]),Normal(),
                IdentityLink(),wts=w2[:])
                dtt[iband] = coef(model)[1]
                err[iband] = stderror(model)[1]
            else
                print("not enough points to estimate dv/v")
                dtt[iband]=NaN
                err[iband]=NaN
            end
        end
    end

    return freqbands, dtt.*100, err.*100
end

"""
Compute wavelet cross spectrum given an array of frequencies
"""
function wxs_allfreq(cur, ref, t, twindow, fwindow, dj, s0, J; wvn::String="morlet", unwrapflag::Bool=false, sig::Bool=false, siglvl::Float64=0.95, windowflag::Bool=true, nwindow::Float64=1.5)
    # Part 1 ----- Cross spectrum analysis
    # define sample frequency
    dt = t[2] - t[1]
    fs = 1/dt
    (fmin, fmax) = fwindow

    # perform wavelet coherence transform
    WXS, WXA, WCT, aWCT, coi, freqs, sig = wct_modified(cur, ref, dt, dj, s0, J, wavelet=wvn, sig=sig, siglvl=siglvl)

    # do inverse cwt for different frequency bands
    if unwrapflag==true
        phase = unwrap(aWCT, dims=2)
    else
        phase = aWCT
    end

    delta_t = Array{Float64,2}(undef, size(phase))
    for i=1:size(phase,2)
        delta_t[:, i] = phase[:, i] ./ (2π.*freqs)
    end

    # zero out cone of influence and data outside frequency band
    if (fmax > maximum(freqs) || fmax<=fmin)
        println("Error: please input correct frequency limits in the frequency window!")
    else
        freq_indout = findall(f->(f < fmin || f > fmax), freqs)
        delta_t[freq_indout, :] .= 0.
        freq_indin = findall(f->(f >= fmin && f <= fmax), freqs)
    end

    # Part 2 ------ dv/v
    dvv = zeros(size(freqs))
    err = zeros(size(freqs))

    if windowflag
        (tmin, tmax) = twindow

        # time checks
        if tmin < minimum(twindow) || tmax > maximum(twindow) || tmax <= tmin
            println("Error: please input right time limits in the time window!")
        else
            # trim time vector
            t_ind = findall(x->(x≤tmax && x≥tmin), t)
            # trim vectors to correspond to time
            wt = t[t_ind]
        end

        # perform linear regression
        for ifreq in freq_indin
            if length(wt) > 2
                if delta_t[ifreq, :] == nothing
                    continue
                end

                # use WXA as weight for regression
                w = 1 ./ WCT[ifreq, :]
                infNaN = findall(x->(isnan.(x) || isinf.(x)), w)
                if length(infNaN)!=0
                    w[infNaN] .= 1.0
                end
                # perform linear regression
                model = glm(@formula(Y ~0 + X),DataFrame(X=wt,Y=delta_t[ifreq,:]),Normal(),
                            IdentityLink(),wts=w)

                dvv[ifreq] = coef(model)[1]
                err[ifreq] = stderror(model)[1]
            else
                println("Not enough points to estimate dv/v")
                dvv[ifreq] = NaN
                err[ifreq] = NaN
            end
        end
    else # using dynamic time window
        for ifreq in freq_indin
            tmin = twindow[1]
            tmax = twindow[1] + nwindow*(1.0 / freqs[ifreq])

            # time checks
            if tmin < minimum(twindow) || tmax > maximum(twindow) || tmax <= tmin
                println("Error: please input right time limits in the time window!")
            else
                # trim time vector
                t_ind = findall(x->(x≤tmax && x≥tmin), t)
                # trim vectors to correspond to time
                wt = t[t_ind]
            end

            if length(wt) > 2
                if delta_t[ifreq, :] == nothing
                    continue
                end

                # use WXA as weight for regression
                w = 1 ./ WCT[ifreq, :]
                infNaN = findall(x->(isnan.(x) || isinf.(x)), w)
                if length(infNaN)!=0
                    w[infNaN] .= 1.0
                end
                # perform linear regression
                model = glm(@formula(Y ~0 + X),DataFrame(X=wt,Y=delta_t[ifreq,:]),Normal(),
                            IdentityLink(),wts=w)

                dvv[ifreq] = coef(model)[1]
                err[ifreq] = stderror(model)[1]
            else
                println("Not enough points to estimate dv/v")
                dvv[ifreq] = NaN
                err[ifreq] = NaN
            end
        end
    end

    return freqs, dvv.*100, err.*100
end

function wts_freqbands(cur, ref, t, twindow, freqbands, dj, s0, J; dvmax::Float64=0.1, ndv::Int64=100, wvn::String="morlet", normalize::Bool=true, windowflag::Bool=true, nwindow::Float64=1.5)
    # define sample frequency
    dt = t[2] - t[1]
    fs = 1/dt

    # apply cwt to two traces
    # define wavelet
    if wvn=="morlet"
        wav = WT.Morlet(6)
    end
    # calculate the CWT of the time series, using identical parameters for both calculations
    cwt1, sj, freqs, coi = cwt(cur, CFW(wav,1/dj), J1=J, dt=dt, s0=s0)
    cwt2, sj, freqs, coi = cwt(ref, CFW(wav,1/dj), J1=J, dt=dt, s0=s0)

    (nband, nfreq) = size(freqbands)
    if nband==0 || nfreq != 2
        println("Error: please check if inputs in freqbands is right!")
    else
        dvv = zeros(Float32, nband)
        cc = zeros(Float32, nband)
        cdp = zeros(Float32, nband)
        err = zeros(Float32, nband)

        # iterate over frequency bands
        for iband=1:nband
            (fmin, fmax) = freqbands[iband, :]

            if (fmax > maximum(freqs) || fmax <= fmin)
                println("Error: please ensure columns 1 and 2 are right frequency limits in freqbands!")
            else
                freq_ind = findall(f->(f>=fmin && f<=fmax), freqs)

                # perform icwt
                icwt1 = icwt(cwt1[freq_ind, :], CFW(wav,1/dj), sj[freq_ind], dt=dt, dj=dj)
                icwt2 = icwt(cwt2[freq_ind, :], CFW(wav,1/dj), sj[freq_ind], dt=dt, dj=dj)

                if windowflag
                    tmin = twindow[1]
                    tmax = twindow[2]

                    # time checks
                    if tmin < minimum(t) || tmax > maximum(t) || tmax <= tmin
                        println("Error: please input right time limits in the time window!")
                    else
                        # trim time vector
                        t_ind = findall(x->(x≤tmax && x≥tmin), t)
                    end
                else
                    tmin = twindow[1]
                    tmax = twindow[1] + nwindow*(1/fmin + 1/fmax)/2.

                    # time checks
                    if tmin < minimum(t) || tmax > maximum(t) || tmax <= tmin
                        println("Error: please input right time limits in the time window!")
                    else
                        # trim time vector
                        t_ind = findall(x->(x≤tmax && x≥tmin), t)
                    end
                end
                # trim vectors to correspond to time
                wt = t[t_ind]
                window = collect(1:length(wt))
                rcwt1 = real.(icwt1)
                rcwt2 = real.(icwt2)
                wcwt1 = rcwt1[t_ind]
                wcwt2 = rcwt2[t_ind]

                # normalize both signals, if appropriate
                if normalize
                    ncwt1 = ((wcwt1 .- mean(wcwt1)) ./ std(wcwt1))[:]
                    ncwt2 = ((wcwt2 .- mean(wcwt2)) ./ std(wcwt2))[:]
                else
                    ncwt1 = wcwt1[:]
                    ncwt2 = wcwt2[:]
                end

                # use stretching to extract dv/v
                (dvv[iband], cc[iband], cdp[iband], eps, err[iband], allC) = stretching(ncwt1, ncwt2, wt, window, fmin, fmax, dvmin=-dvmax, dvmax=dvmax)
            end
        end
    end

    return freqbands, dvv, cc, cdp, err
end

function wts_allfreqs(cur, ref, t, twindow, fwindow, dj, s0, J; dvmax::Float64=0.1, ndv::Int64=100, wvn::String="morlet", normalize::Bool=true, windowflag::Bool=true, nwindow::Float64=1.5)
    # define sample frequency
    dt = t[2] - t[1]
    fs = 1/dt
    (fmin, fmax) = fwindow

    # apply cwt to two traces
    # define wavelet
    if wvn=="morlet"
        wav = WT.Morlet(6)
    end
    # calculate the CWT of the time series, using identical parameters for both calculations
    cwt1, sj, freqs, coi = cwt(cur, CFW(wav,1/dj), J1=J, dt=dt, s0=s0)
    cwt2, sj, freqs, coi = cwt(ref, CFW(wav,1/dj), J1=J, dt=dt, s0=s0)

    # extract real part of cwt
    rcwt1 = real.(cwt1)
    rcwt2 = real.(cwt2)

    # zero out cone of influence and data outside frequency band
    if (fmax > maximum(freqs) || fmax<=fmin)
        println("Error: please input correct frequency limits in the frequency window!")
    else
        freq_indout = findall(f->(f < fmin || f > fmax), freqs)
        rcwt1[freq_indout, :] .= 0.
        rcwt2[freq_indout, :] .= 0.
        freq_indin = findall(f->(f >= fmin && f <= fmax), freqs)

        # Use stretching method to extract dvv
        nfreq=length(freqs)

        dvv = zeros(Float32, nfreq)
        cc = zeros(Float32, nfreq)
        cdp = zeros(Float32, nfreq)
        err = zeros(Float32, nfreq)

        for ifreq in freq_indin
            if windowflag
                tmin = twindow[1]
                tmax = twindow[2]

                # time checks
                if tmin < minimum(t) || tmax > maximum(t) || tmax <= tmin
                    println("Error: please input right time limits in the time window!")
                else
                    # trim time vector
                    t_ind = findall(x->(x≤tmax && x≥tmin), t)
                end
            else
                tmin = twindow[1]
                tmax = twindow[1] + nwindow*(1/fmin + 1/fmax)/2.

                # time checks
                if tmin < minimum(t) || tmax > maximum(t) || tmax <= tmin
                    println("Error: please input right time limits in the time window!")
                else
                    # trim time vector
                    t_ind = findall(x->(x≤tmax && x≥tmin), t)
                end
            end

            # prepare time axis and its indices
            wt = t[t_ind]
            it = collect(1:length(wt))
            # prepare windowed data
            wcwt1 = rcwt1[ifreq, t_ind]
            wcwt2 = rcwt2[ifreq, t_ind]

            # normalize both signals if appropriate
            # normalize both signals, if appropriate
            if normalize
                ncwt1 = ((wcwt1 .- mean(wcwt1)) ./ std(wcwt1))[:]
                ncwt2 = ((wcwt2 .- mean(wcwt2)) ./ std(wcwt2))[:]
            else
                ncwt1 = wcwt1[:]
                ncwt2 = wcwt2[:]
            end

            # compute the stretching
            window = collect(1:length(wt))
            # use stretching to extract dv/v
            (dvv[ifreq], cc[ifreq], cdp[ifreq], eps, err[ifreq], allC) = stretching(ncwt1, ncwt2, wt, window, fmin, fmax, dvmin=-dvmax, dvmax=dvmax)
        end
    end

    return freqs, dvv, cc, cdp, err
end

function wtdtw_freqbands(cur, ref, t, twindow, freqbands, maxLag, b, direction, dj, s0, J; wvn::String="morlet", normalize::Bool=true, windowflag::Bool=true, nwindow::Float64=1.5)
    # define sample frequency
    dt = t[2] - t[1]
    fs = 1/dt

    # apply cwt to two traces
    # define wavelet
    if wvn=="morlet"
        wav = WT.Morlet(6)
    end
    # calculate the CWT of the time series, using identical parameters for both calculations
    cwt1, sj, freqs, coi = cwt(cur, CFW(wav,1/dj), J1=J, dt=dt, s0=s0)
    cwt2, sj, freqs, coi = cwt(ref, CFW(wav,1/dj), J1=J, dt=dt, s0=s0)

    (nband, nfreq) = size(freqbands)
    if nband==0 || nfreq != 2
        println("Error: please check if inputs in freqbands is right!")
    else
        dvv = zeros(Float32, nband)
        err = zeros(Float32, nband)

        # iterate over frequency bands
        for iband=1:nband
            (fmin, fmax) = freqbands[iband, :]

            if (fmax > maximum(freqs) || fmax <= fmin)
                println("Error: please ensure columns 1 and 2 are right frequency limits in freqbands!")
            else
                freq_ind = findall(f->(f>=fmin && f<=fmax), freqs)

                # perform icwt
                icwt1 = icwt(cwt1[freq_ind, :], CFW(wav,1/dj), sj[freq_ind], dt=dt, dj=dj)
                icwt2 = icwt(cwt2[freq_ind, :], CFW(wav,1/dj), sj[freq_ind], dt=dt, dj=dj)

                if windowflag
                    tmin = twindow[1]
                    tmax = twindow[2]

                    # time checks
                    if tmin < minimum(t) || tmax > maximum(t) || tmax <= tmin
                        println("Error: please input right time limits in the time window!")
                    else
                        # trim time vector
                        t_ind = findall(x->(x≤tmax && x≥tmin), t)
                    end
                else
                    tmin = twindow[1]
                    tmax = twindow[1] + nwindow*(1/fmin + 1/fmax)/2.

                    # time checks
                    if tmin < minimum(t) || tmax > maximum(t) || tmax <= tmin
                        println("Error: please input right time limits in the time window!")
                    else
                        # trim time vector
                        t_ind = findall(x->(x≤tmax && x≥tmin), t)
                    end
                end
                # trim vectors to correspond to time
                wt = t[t_ind]
                window = collect(1:length(wt))
                rcwt1 = real.(icwt1)
                rcwt2 = real.(icwt2)
                wcwt1 = rcwt1[t_ind]
                wcwt2 = rcwt2[t_ind]

                # normalize both signals, if appropriate
                if normalize
                    ncwt1 = ((wcwt1 .- mean(wcwt1)) ./ std(wcwt1))[:]
                    ncwt2 = ((wcwt2 .- mean(wcwt2)) ./ std(wcwt2))[:]
                else
                    ncwt1 = wcwt1[:]
                    ncwt2 = wcwt2[:]
                end

                (stbarTime, stbar, dist, error) = dtwdt(ncwt1, ncwt2, dt, maxLag=maxLag, b=b, direction=direction)

                # perform linear regression
                model = glm(@formula(Y ~0 + X),DataFrame(X=wt,Y=stbarTime),Normal(),
                            IdentityLink(),wts=ones(length(wt)))

                dvv[iband] = coef(model)[1]
                err[iband] = stderror(model)[1]
            end
        end
    end
    return freqs, dvv, err
end

function wtdtw_allfreqs(cur, ref, t, twindow, fwindow, maxLag, b, direction, dj, s0, J; wvn::String="morlet", normalize::Bool=true, windowflag::Bool=true, nwindow::Float64=1.5)
    # define sample frequency
    dt = t[2] - t[1]
    fs = 1/dt
    (fmin, fmax) = fwindow

    # apply cwt to two traces
    # define wavelet
    if wvn=="morlet"
        wav = WT.Morlet(6)
    end
    # calculate the CWT of the time series, using identical parameters for both calculations
    cwt1, sj, freqs, coi = cwt(cur, CFW(wav,1/dj), J1=J, dt=dt, s0=s0)
    cwt2, sj, freqs, coi = cwt(ref, CFW(wav,1/dj), J1=J, dt=dt, s0=s0)

    # extract real part of cwt
    rcwt1 = real.(cwt1)
    rcwt2 = real.(cwt2)

    # zero out cone of influence and data outside frequency band
    if (fmax > maximum(freqs) || fmax<=fmin)
        println("Error: please input correct frequency limits in the frequency window!")
    else
        freq_indout = findall(f->(f < fmin || f > fmax), freqs)
        rcwt1[freq_indout, :] .= 0.
        rcwt2[freq_indout, :] .= 0.
        freq_indin = findall(f->(f >= fmin && f <= fmax), freqs)

        # Use stretching method to extract dvv
        nfreq=length(freqs)

        dvv = zeros(Float32, nfreq)
        err = zeros(Float32, nfreq)

        for ifreq in freq_indin
            if windowflag
                tmin = twindow[1]
                tmax = twindow[2]

                # time checks
                if tmin < minimum(t) || tmax > maximum(t) || tmax <= tmin
                    println("Error: please input right time limits in the time window!")
                else
                    # trim time vector
                    t_ind = findall(x->(x≤tmax && x≥tmin), t)
                end
            else
                tmin = twindow[1]
                tmax = twindow[1] + nwindow*(1/fmin + 1/fmax)/2.

                # time checks
                if tmin < minimum(t) || tmax > maximum(t) || tmax <= tmin
                    println("Error: please input right time limits in the time window!")
                else
                    # trim time vector
                    t_ind = findall(x->(x≤tmax && x≥tmin), t)
                end
            end

            # prepare time axis and its indices
            wt = t[t_ind]
            it = collect(1:length(wt))
            # prepare windowed data
            wcwt1 = rcwt1[ifreq, t_ind]
            wcwt2 = rcwt2[ifreq, t_ind]

            # normalize both signals if appropriate
            # normalize both signals, if appropriate
            if normalize
                ncwt1 = ((wcwt1 .- mean(wcwt1)) ./ std(wcwt1))[:]
                ncwt2 = ((wcwt2 .- mean(wcwt2)) ./ std(wcwt2))[:]
            else
                ncwt1 = wcwt1[:]
                ncwt2 = wcwt2[:]
            end

            (stbarTime, stbar, dist, error) = dtwdt(ncwt1, ncwt2, dt, maxLag=maxLag, b=b, direction=direction)

            # perform linear regression
            model = glm(@formula(Y ~0 + X),DataFrame(X=wt,Y=stbarTime),Normal(),
                        IdentityLink(),wts=ones(length(wt)))

            dvv[ifreq] = coef(model)[1]
            err[ifreq] = stderror(model)[1]
        end
    end

    return freqs, dvv, err
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
abstract type AbstractElementGrid{T} <: AbstractVector{T} end
abstract type AbstractBasis{T} <: AbstractVector{T} end

struct Element1D
    a::Float64
    b::Float64
    function Element1D(a,b)
        @argcheck a < b
        new(a,b)
    end
end

struct ElementGrid{T} <: AbstractElementGrid{T}
    el::Element1D
    gl::GaussLegendre{T}
    scaling::T
    shift::T
    function ElementGrid(a, b, n)
        new{Float64}(Element1D(a,b), GaussLegendre(n), (b-a)/2, (a+b)/2)
    end
end

struct ElementGridLobatto{T} <: AbstractElementGrid{T}
    el::Element1D
    gl::LobattoLegendre{T}
    scaling::T
    shift::T
    function ElementGridLobatto(a, b, n)
        new{Float64}(Element1D(a,b), LobattoLegendre(n), (b-a)/2, (a+b)/2)
    end
end

struct Basis{T} <: AbstractBasis{T}
    egvector::Vector{ElementGrid{T}}
    x::Vector{T}
    function Basis(eg::ElementGrid{T}...) where T
        if length(eg) > 1
            for i in 2:length(eg)
                @assert eg[i-1] <= eg[i]
            end
        end
        new{T}(collect(eg), vcat(eg...))
    end
end

struct BasisLobatto{T} <: AbstractBasis{T}
    egvector::Vector{ElementGridLobatto{T}}
    x::Vector{T}
    D::Matrix{T}
    w::Vector{T}
    function BasisLobatto(eg::ElementGridLobatto{T}...) where T
        if length(eg) > 1
            x = Vector(eg[1])
            for i in 2:length(eg)
                @assert eg[i-1][end] ≈ eg[i][begin]
                append!(x, eg[i][begin+1:end])
            end

            # Build derivative matrix and integration weights
            dv = [ derivative_matrix(x) for x in eg ]
            cr = [size(x,1) for x in dv] # dv is square matrix
            # Size of derivative matrix. Lobatto basis overlaps with adjacent elements
            # and we need to remove this overlap, by reducing number of points.
            nd = sum(cr) - ( length(cr) - 1 )
            D = zeros(T ,nd, nd)  # Derivate matrix
            w = zeros(T, nd)      # Integration weights
            D[1:cr[1], 1:cr[1]] = dv[1]
            w[1:cr[1]] = get_weight(eg[1])
            # Index range for final derivative matrix build from element ones
            tmp = ones(Int, length(cr))
            tmp[1] = 0
            cum_index = cumsum( cr .- tmp )
            for i in 2:length(cum_index)
                ir = cum_index[i-1]:cum_index[i]
                D[ir, ir] += dv[i]
                D[cum_index[i-1],:] .*= 0.5 # remove double count from overlap point
                w[ir] += get_weight(eg[i])
            end
        else
            x = Vector(eg[1])
            D = eg[1].gl.D
            w = get_weight(eg[1])
        end
        new{T}(collect(eg), x, D, w)
    end
end

function Basis(a, b, ne, np)
    @argcheck a < b
    step = (b-a)/ne
    eg = [ElementGrid( a+(i-1)*step, a + i*step, np )  for i in 1:ne]
    return Basis(eg...)
end

function BasisLobatto(a, b, ne, np)
    @argcheck a < b
    step = (b-a)/ne
    eg = [ElementGridLobatto( a+(i-1)*step, a + i*step, np )  for i in 1:ne]
    return BasisLobatto(eg...)
end

function Base.show(io::IO, ::MIME"text/plain", eg::ElementGrid)
    print(io, "ElementGrid - $(length(eg)) points")
end

function Base.show(io::IO, ::MIME"text/plain", b::Basis)
    print(io, "Basis - $(length(b.egvector)) elements and $(length(b)) points")
end

Base.size(eg::AbstractElementGrid) = size(eg.gl.nodes)
Base.size(b::AbstractBasis) = size(b.x)

Base.firstindex(b::AbstractBasis) = firstindex(b.x)
Base.firstindex(eg::AbstractElementGrid) = firstindex(eg.gl.nodes)

Base.lastindex(b::AbstractBasis) = lastindex(b.x)
Base.lastindex(eg::AbstractElementGrid) = lastindex(eg.gl.nodes)

Base.getindex(eg::AbstractElementGrid, i::Int) = eg.gl.nodes[i]*eg.scaling + eg.shift
Base.getindex(b::AbstractBasis, i::Int) = b.x[i]

for OP in [:(<), :(<=)]
    @eval begin
        function Base.$OP(a::Element1D, b::Element1D)
            return $OP(a.b, b.a)
        end
        function Base.$OP(a::AbstractElementGrid, b::AbstractElementGrid)
            return $OP(a.el, b.el)
        end
    end
end

for OP in [:(>), :(>=)]
    @eval begin
        function Base.$OP(a::Element1D, b::Element1D)
            return $OP(a.a, b.b)
        end
        function Base.$OP(a::AbstractElementGrid, b::AbstractElementGrid)
            return $OP(a.el, b.el)
        end
    end
end

function get_weight(eg::AbstractElementGrid, i::Int)
    return eg.gl.weights[i] * eg.scaling
end

function get_weight(eg::AbstractElementGrid)
    return eg.gl.weights .* eg.scaling
end

function PolynomialBases.derivative_matrix(eg::AbstractElementGrid)
    return eg.gl.D ./ eg.scaling
end

function metric_tensor(b::AbstractBasis)
    return Diagonal(get_weight(b))
end

function PolynomialBases.derivative_matrix(b::Basis)
    dv = [ derivative_matrix(x) for x in b.egvector  ]
    # get block sizes for block banded matrix
    cr = [size(x)[1] for x in dv] # dv is square matrix
    out = BlockBandedMatrix{eltype(b)}(undef, cr, cr, (0,0))
    for i in 1:length(cr)
        out[Block(i,i)] = dv[i]
    end
    return Matrix(out)
end

function PolynomialBases.derivative_matrix(b::BasisLobatto)
    return b.D
end

"""
    get_weight(b::AbstractBasis)

Integration weights for the basis
"""
function get_weight(b::Basis)
    tmp = get_weight.(b.egvector)
    return vcat(tmp...)
end

function get_weight(b::BasisLobatto)
    return b.w
end

"""
    get_identity(b::AbstractBasis)

Identity matrix for the given basis
"""
function get_identity(b::AbstractBasis)
    return diagm(ones(length(b)))
end

"""
    get_length(el::Element1D)
    get_length(eg::ElementGrid)
    get_length(b::Basis)

Gives physical length the input represents
"""
function get_length(el::Element1D)
    return el.b - el.a
end

function get_length(eg::AbstractElementGrid)
    return get_length(eg.el)
end

function get_length(b::AbstractBasis)
    return sum( get_length, b.egvector)
end

function get_element(b::AbstractBasis, i)
    return b.egvector[i]
end

function number_of_elements(b::AbstractBasis)
    return length(b.egvector)
end

function get_element_indexes(b::BasisLobatto, i::Int)
    @argcheck 0 < i <= number_of_elements(b)
    l1 = sum( [ length( get_element(b, x) ) for x in 1:i ] ) - max(0, i-1)
    l0 = sum( [ length( get_element(b, x) ) for x in 1:i-1 ] ) - max(0, i-1) + 1
    return l0:l1
end

function get_element_indexes(b::BasisLobatto, i)
    return (sort!∘union)( [ get_element_indexes(b, x) for x in i ]...  )
end

max_stars_repo_name: tjjarvinen/OneDimensionalCorrelation

max_issues_repo_name: tjjarvinen/OneDimensionalCorrelation

max_forks_repo_name: tjjarvinen/OneDimensionalCorrelation


text:
using SIMDHints: ivdep
using BenchmarkTools

f!(y, x) = y .= ivdep.(x .+ y)
g!(y, x) = y .= x .+ y

a = randn(2^10)
y = zero(a)

@btime f!($y, $y) setup=(fill!($y, 0))
@btime g!($y, $y) setup=(fill!($y, 0))
nothing


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_ps_opt_rho1_9(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_ps_opt_rho1_9!(A_copy)
end

@inline function exp_ps_opt_rho1_9!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=8
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 T2k8
    # Computing Bb2 = x*I+x*A
    coeff1=0.30083888547588206
    coeff2=1.0751360702391695
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.30083888547588206
    coeff2=1.0751360702391695
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.3287994669675759
    coeff2=0.6302693938160198
    coeff3=1.536487303722
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=0.6291288319953808
    coeff2=1.5445527760735984
    coeff3=0.06828138863920057
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=0.01820773746059349
    coeff2=0.36143000551624754
    coeff3=0.7785296989138268
    coeff4=1.0178724684635254
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=-3.901563650729534
    coeff2=0.29276483269826203
    coeff3=0.02397633499952574
    coeff4=0.0008026495355626526
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=6.97556294530088e-7
    coeff2=0.011883016278515129
    coeff3=0.3610383308366146
    coeff4=-4.819055642175235
    coeff5=0.6031240570850653
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-4.763068434972627e-6
    coeff2=1.2563376565457575e-7
    coeff3=9.983587156878525e-9
    coeff4=5.701593011031481e-10
    coeff5=1.527146310377267e-11
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=-0.003910315472289407
    coeff2=-0.04148590496173931
    coeff3=0.8367416720923179
    coeff4=1.2101133147668501
    coeff5=0.9887901240572256
    coeff6=0.003241259889641698
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.22742373870204555
    coeff2=0.00441487179797981
    coeff3=0.0005139372684020878
    coeff4=3.036579116890275e-5
    coeff5=1.0042397964978794e-5
    coeff6=0.9999999134512341
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing T2k8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.9682967326425941
    coeff2=0.8510597305281336
    coeff3=0.3004272311247782
    coeff4=0.039414194971706944
    coeff5=-0.15919544727149795
    coeff6=-4.2294712847109625e-6
    coeff7=1.0207974923352356
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    return memslots1 # Returning T2k8
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -0.7481803354345802 <= q <= 3.1247048530266452)

                     add_NL_constraint(m, :(softplus(0.30283339323429637 + -0.6237495659853991*softplus(-0.21938578694289834 + -0.903262746252536*softplus(0.9051056598158431 + -0.13053580855908242*$(x[1]) + -0.5492925144525009*$(x[2]) + 0.019096353213378325*$(x[3]) + 0.3378807953855012*$(x[4])) + 0.16468819845671145*softplus(0.06066012774712304 + 0.010865522181132903*$(x[1]) + -0.27940603676198705*$(x[2]) + -0.359831046364532*$(x[3]) + -0.7556261576696839*$(x[4]))) + 0.9377177734653888*softplus(0.5819650565710264 + -0.47161154460717736*softplus(0.9051056598158431 + -0.13053580855908242*$(x[1]) + -0.5492925144525009*$(x[2]) + 0.019096353213378325*$(x[3]) + 0.3378807953855012*$(x[4])) + -0.23723822344605994*softplus(0.06066012774712304 + 0.010865522181132903*$(x[1]) + -0.27940603676198705*$(x[2]) + -0.359831046364532*$(x[3]) + -0.7556261576696839*$(x[4])))) + softplus(-0.5162811028857681 + 0.4410139812710443*softplus(-0.21938578694289834 + -0.903262746252536*softplus(0.9051056598158431 + -0.13053580855908242*$(x[1]) + -0.5492925144525009*$(x[2]) + 0.019096353213378325*$(x[3]) + 0.3378807953855012*$(x[4])) + 0.16468819845671145*softplus(0.06066012774712304 + 0.010865522181132903*$(x[1]) + -0.27940603676198705*$(x[2]) + -0.359831046364532*$(x[3]) + -0.7556261576696839*$(x[4]))) + 0.2178047619382908*softplus(0.5819650565710264 + -0.47161154460717736*softplus(0.9051056598158431 + -0.13053580855908242*$(x[1]) + -0.5492925144525009*$(x[2]) + 0.019096353213378325*$(x[3]) + 0.3378807953855012*$(x[4])) + -0.23723822344605994*softplus(0.06066012774712304 + 0.010865522181132903*$(x[1]) + -0.27940603676198705*$(x[2]) + -0.359831046364532*$(x[3]) + -0.7556261576696839*$(x[4])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#function rrule(::typeof(MPO), x::Vector{<:ITensor}; kwargs...)
#  y = MPO(x; kwargs...)
#  #@show y
#  function MPO_pullback(ȳ)
#    #@show ȳ
#    return ȳ.data
#    #ȳtensors = ȳ.data
#    #n = length(ȳtensors)
#    #envL = [ȳtensors[1] * dag(x[1]), ]
#    #envR = [ȳtensors[n] * dag(x[n]), ]
#    #for j in 2:n-1
#    #  push!(envL, envL[j-1] * ȳtensors[j] * dag(x[j]))
#    #  push!(envR, envR[j-1] * ȳtensors[n+1-j] * dag(x[n+1-j]))
#    #end
#    #x̄= ITensor[]
#    #push!(x̄, ȳtensors[1]  * envR[n-1])
#    #for j in 2:n-1
#    # push!(x̄, envL[j-1] * ȳtensors[j] * envR[n-j])
#    #end
#    #push!(x̄, envL[n-1] * ȳtensors[n])
#    #return (NoTangent(), x̄)
#  end
#  return y, MPO_pullback
#end

function rrule(::typeof(*), x1::MPO, x2::MPO; kwargs...)
  y = *(x1, x2; kwargs...)
  function contract_pullback(ȳ)
    x̄1 = *(ȳ, dag(x2); kwargs...)
    x̄2 = *(dag(x1), ȳ; kwargs...)
    return (NoTangent(), x̄1, x̄2)
  end
  return y, contract_pullback
end

function rrule(::typeof(tr), x::MPO; kwargs...)
  y = tr(x; kwargs...)
  function contract_pullback(ȳ)
    s = noprime(firstsiteinds(x))
    n = length(s)
    x̄ = ȳ * MPO(s, "Id")

    plev = get(kwargs, :plev, 0 => 1)
    for j in 1:n
      x̄[j] = mapprime(x̄[j], 0 => first(plev), 1 => last(plev))
    end
    return (NoTangent(), ȳ * x̄)
  end
  return y, contract_pullback
end

function rrule(::typeof(inner), x1::MPS, x2::MPO, x3::MPS; kwargs...)
  if !hassameinds(siteinds, x1, (x2, x3)) || !hassameinds(siteinds, x3, (x2, x1))
    error(
      "Taking gradients of `inner(x::MPS, A::MPO, y::MPS)` is not supported if the site indices of the input MPS and MPO don't match. Try using if you input `inner(x, A, y), try `inner(x', A, y)` instead.",
    )
  end

  y = inner(x1, x2, x3; kwargs...)
  function inner_pullback(ȳ)
    x̄1 = dag(ȳ) * contract(x2, x3; kwargs...)
    x̄2 = ȳ * dag(_contract(MPO, dag(x1), x3; kwargs...))
    x̄3 = contract(dag(x2), x1; kwargs...) * ȳ

    @assert siteinds(x1) == siteinds(x̄1)
    @assert hassameinds(siteinds, x2, x̄2)
    @assert siteinds(x3) == siteinds(x̄3)

    return (NoTangent(), x̄1, x̄2, x̄3)
  end
  return y, inner_pullback
end


text:
using LinearAlgebra
using Statistics
using DelimitedFiles;
using jInv.Mesh;
using MAT;
using jInv
using jInv.InverseSolve
using ShapeReconstructionPaLS
using ShapeReconstructionPaLS.DipTransform;
using ParamLevelSet;
using ShapeReconstructionPaLS.ShapeFromSilhouette;
using ShapeReconstructionPaLS.Utils;
using Statistics
using Distributed
using SparseArrays
using Random

ENV["MPLBACKEND"] = "Qt4Agg"

plotting = false;

if plotting
	using PyPlot;
	using jInvVis;
	close("all");
end

####
#Parameters to choose if: inversion with dip only/visual hull only/ joint inversion
invertDip = true;
invertVis = false;
#####

#Place new added RBFs according to gradient values or randomly:
locateRBFwithGrads = true;

################################################
### Reading the model and pad ##################
################################################
model = "skull_160"; ## skull
 #model = "dancing220";
file = matopen(string(pwd(),"/models/",model,".mat"));

A = read(file, "B");
close(file);
A = convert(Array{Float32,3},A[1:1:end,1:1:end,1:1:end]);
n_inner = collect(size(A));
pad = ceil.(Int64,0.1*n_inner[1]);
n = n_inner .+ 2*pad;

n_tup = tuple(n...)
volume = prod(n);
n1 = n_tup[1];
n2 = n_tup[2];
n3 = n_tup[3];

m = zeros(Float64,n_tup);
mask = zeros(Int8,n_tup);
m[pad+1:end-pad,pad+1:end-pad,pad+1:end-pad] = A; A = 0;
mask[pad+1:end-pad,pad+1:end-pad,pad+1:end-pad] .= 1;
if plotting
	plotModel(m);
	savefig("TrueModel.png");
end

################################################################
### Preparing the synthetic data using mat-free implementation
################################################################
Mesh = getRegularMesh([0.0 5.0 0.0 5.0 0.0 5.0],n);
ScreenMesh = getRegularMesh(Mesh.domain[3:6],n[2:3]); ## for vis.
ScreenDomain = ScreenMesh.domain;
n_screen = ScreenMesh.n;

noiseAnglesDeg =0.0#4.0;
#original was mean(2*Mesh.h);
noiseTrans =0.0#mean(8*Mesh.h);
noiseAngles =deg2rad(noiseAnglesDeg);
#########################################################################################################
### Dipping Data ########################################################################################
#########################################################################################################
samplingBinning = 2;
nDips = 0;
nShots = 0;
if invertDip
	println("Data gen resolution:")
	println(n)
	noiseSample = 2.0*prod(Mesh.h);
	nDips = 100;
	theta_phi_dip = deg2rad.(convert(Array{Float64,2}, [rand(0:359,nDips)  rand(0:90,nDips)]));
	dipDataFilename = string("dipData_",model,"_dataRes",n,"_bin",samplingBinning,"_noiseSampleAnglesTrans",[noiseSample;noiseAngles;noiseTrans],"_dips",nDips);
	prepareSyntheticDipDataFiles(m,Mesh,dipDataFilename,pad,theta_phi_dip,samplingBinning, noiseSample, noiseAngles,noiseTrans,MATFree);

	#########################################################################################################
	### Read the dipping data ###############################################################################
	#########################################################################################################
	n,domain,samplingBinning,DataDip,theta_phi_dip,pad = readDipDataFile(dipDataFilename);
end

####################################

#################################################################################################################
### Generate and Read the visual data ###############################################################################
#################################################################################################################

if invertVis
	nShots = 16;
	noiseSample = 0.0;
	factorShotsCreation = 4;
	theta_phi_vis = deg2rad.(convert(Array{Float64,2}, [rand(0:359,factorShotsCreation*nShots)  rand(0:90,factorShotsCreation*nShots)]));

	VisDataFilename = string("visData_",model,"_dataRes",n,"_noiseSampleAnglesTrans",[noiseSample;noiseAngles;noiseTrans],"_shots",nShots);

	XlocScreen = Mesh.domain[2];
	LocCamera = [-16.0*(Mesh.domain[2] - Mesh.domain[1]) ; 0.0 ;0.0]; ## 0.0 here is the middle of the domain
	noiseCamLoc = 0.0;
	noiseScreenLoc = 0.0;
	prepareSyntheticVisDataFiles(m,Mesh,ScreenMesh,pad,theta_phi_vis,XlocScreen,LocCamera,
	VisDataFilename,factorShotsCreation*nShots, noiseCamLoc, noiseScreenLoc, noiseAngles,noiseTrans,MATFree);

	##################################################################################################################
	##################################################################################################################

	n,domain,XScreenLoc,CameraLoc,ScreenDomain,n_screen,DataVis,theta_phi_vis = readVisDataFile(VisDataFilename);
	norms = vec(sum(DataVis,1));
	idsorted = sortperm(norms);
	theta_phi_vis = theta_phi_vis[idsorted[(end-nShots+1):end],:];
	DataVis = DataVis[:,idsorted[(end-nShots+1):end]];
	# for kkk = 1:3:nShots
		# dk = DataVis[:,kkk];
		# dk = reshape(dk,tuple(n_screen...));
		# figure();
		# plotModel(dk)
	# end

end
ScreenMesh 	= getRegularMesh(ScreenDomain,div.(n_screen,samplingBinning));
Mesh = getRegularMesh(domain,div.(n,samplingBinning));
pad = div(pad,samplingBinning);
if samplingBinning==2 && invertVis
	DataVis = coarsenByTwo(DataVis,ScreenMesh.n,nShots);
	if plotting
		figure();
		for kkk = 1:max(4,nShots)
			dk = DataVis[:,kkk];
			dk = reshape(dk,tuple(ScreenMesh.n...));
			subplot(8,8,kkk);
			plotModel(dk);

		end
		savefig(string("VisTrueModel.png"));
	end
end

samplingBinning = 1;
println("Inv mesh:")
println(Mesh.n)

# method = RBFBasedSimple2;
# method = RBF10BasedSimple2;
# method = RBF10BasedSimple1
# method = RBFBasedSimple1

method = RBF10Based; methodName="RBF10Based";

#method = MATFree; methodName="MATFree"

n_tup = tuple(Mesh.n...);
nWorkers = nworkers();
misfun = SSDFun; ## least squares

if invertDip
	### Set up the dip inversion
	pForDip = getDipParam(Mesh,theta_phi_dip,zeros(nDips,3),samplingBinning,nWorkers,method);

	### Create Dip param if we use the GRAD of MATFREE to locate new RBFs
	if(locateRBFwithGrads)
		pForDip_MATFree = getDipParam(Mesh,theta_phi_dip,zeros(nDips,3),samplingBinning,nWorkers,MATFree);

	end
	dobsDip = divideDataToWorkersData(nWorkers,DataDip);
	Wd_Dip = Array{Array{Float32}}(undef,nWorkers);
	for k=1:length(Wd_Dip)
		Wd_Dip[k] = ones(Float32,size(dobsDip[k]));
	end

end

if invertVis
### Set up the vis inversion
	pForVis = getVisHullParam(Mesh,theta_phi_vis,zeros(nShots,3),ScreenMesh,XlocScreen,LocCamera,method,nWorkers);
	dobsVis = divideDataToWorkersData(nWorkers,DataVis);
	Wd_Vis = Array{Array{Float32}}(nWorkers);
	for k=1:length(Wd_Vis)
		println("sqrt:",0.1*(Mesh.domain[2]-Mesh.domain[1])*prod(Mesh.h[2:3]));
		Wd_Vis[k] = sqrt(0.05*(Mesh.domain[2]-Mesh.domain[1])*prod(Mesh.h[2:3]))*ones(Float32,size(dobsVis[k]));
		# Wd_Vis[k] = sqrt(0.3)*ones(Float32,size(dobsVis[k]));

	end
end

if invertVis && invertDip
	n_Moves_all = nShots + nDips;
elseif invertDip
	n_Moves_all = nDips;
else
	n_Moves_all = nShots;
end

isSimple = 0;
isRBF10 = 0;

#########################################################################################################
############### Volumetric inversion (pixel-wise) #######################################################
#########################################################################################################
if method == MATBased || method == MATFree
# #### Set active cells
	mask = zeros(Int8,n_tup);
	mask[pad+1:end-pad,pad+1:end-pad,pad+1:end-pad] .= 1;
	Iact = speye(Float32,prod(Mesh.n));
	Iact = Iact[:,mask[:] .== 1];mask = 0;
	IactPlot = Iact;
	sback = zero(Float32);

	#### USE A BOUND MODEL or bounds
	a = 0.0;
	b = 1.0;

	nact = size(Iact,2);
	modfun = identityMod;
	boundsHigh = zeros(Float32,nact) + b;
	boundsLow = zeros(Float32,nact) + a;
	mref = zeros(Float64,nact);
	modfunForPlotting = modfun;

	#### Use smoothness regularization
	regfun = (m, mref, M) -> wTVReg(m, mref, M,Iact=Iact);
	HesPrec=getSSORCGRegularizationPreconditioner(1.0,1e-3,100)
	# HesPrec = getEmptyRegularizationPreconditioner();
	alpha = 1e-10; #1e-10;
	m0 = copy(mref);
###############################################################################################
############## USE RBF dictionary #############################################################
###############################################################################################

else
	nRBF = 5;

	isRBF10 = (method == RBF10BasedSimple1 || method==RBF10BasedSimple2 || method==RBF10Based);
	isSimple = (method==RBFBasedSimple1 || method==RBF10BasedSimple1 || method==RBFBasedSimple2 || method==RBF10BasedSimple2);

	if isRBF10
		numParamOfRBF = 10;
	else
		numParamOfRBF = 5;
	end
	n_m_simple = numParamOfRBF*nRBF;
	nAll = isSimple ? numParamOfRBF*nRBF : numParamOfRBF*nRBF + n_Moves_all*5;

	modfun = identityMod;

	if isSimple
		modfunForPlotting = (m)->ParamLevelSetModelFunc(Mesh,m;computeJacobian = 0,
			sigma = getDefaultHeavySide() ,bf = 1, numParamOfRBF = numParamOfRBF);
	else
		modfunForPlotting = (m)->ParamLevelSetModelFunc(Mesh,m[1:(numParamOfRBF*nRBF)];computeJacobian = 0,
			sigma = getDefaultHeavySide() ,bf = 1, numParamOfRBF = numParamOfRBF);
	end

	Iact = 1.0;
	IactPlot = 1.0;
	sback = zero(Float32);
	a = -100000000.0;
	b = 10000000.0;
	boundsHigh = zeros(Float32,nAll) .+ b;
	boundsLow = zeros(Float32,nAll) .+ a;
	# m0 - initial guess. mref - reference for regularization.
	mref = zeros(Float32,nAll);
	m0   = zeros(Float64,n_m_simple);
	mCenters = zeros(Float32,n_tup);
	mref[1:numParamOfRBF:n_m_simple] = ParamLevelSet.centerHeavySide*2.0*rand(nRBF);
	if isRBF10
		# here we initialize B as an identity matrix
		mref[2:10:n_m_simple] .= 1.0;
		mref[5:10:n_m_simple] .= 1.0;
		mref[7:10:n_m_simple] .= 1.0;
		boundsLow[2:10:n_m_simple] .= 0.05;
		boundsLow[5:10:n_m_simple] .= 0.05;
		boundsLow[7:10:n_m_simple] .= 0.05;
	else
		mref[2:5:n_m_simple] .= 1.0;
		boundsLow[2:5:n_m_simple] .= 0.05;
	end
	mref[(numParamOfRBF-2):numParamOfRBF:n_m_simple] .= (Mesh.domain[1] + Mesh.domain[2])/2;
	mref[(numParamOfRBF-1):numParamOfRBF:n_m_simple] .= (Mesh.domain[3] + Mesh.domain[4])/2;
	mref[numParamOfRBF:numParamOfRBF:n_m_simple] .= (Mesh.domain[5] + Mesh.domain[6])/2;
	m0[:] = mref[1:n_m_simple];
	if !isSimple
		if invertDip
			m0 = wrapRBFparamAndRotationsTranslations(m0,theta_phi_dip,zeros(size(theta_phi_dip,1),3));
		end
		if invertVis
			m0 = wrapRBFparamAndRotationsTranslations(m0,theta_phi_vis,zeros(size(theta_phi_vis,1),3));
		end
	end
	mref[:] = m0;
	m0[(numParamOfRBF-2):numParamOfRBF:n_m_simple] .+= 0.1*(Mesh.domain[2] - Mesh.domain[1])*randn(nRBF);
	m0[(numParamOfRBF-1):numParamOfRBF:n_m_simple] .+= 0.1*(Mesh.domain[4] - Mesh.domain[3])*randn(nRBF);
	m0[numParamOfRBF:numParamOfRBF:n_m_simple]     .+= 0.1*(Mesh.domain[6] - Mesh.domain[5])*randn(nRBF);

	alpha = 5e-1;

	II = (sparse(1.0I, nAll,nAll));
	regfun = (m, mref, M)->TikhonovReg(m,mref,M,II);
	HesPrec = getEmptyRegularizationPreconditioner();

	if isRBF10
		spd_reg = (m,mref,M) -> RBF_SPD_regularization(m,mref,nRBF);
		regfun 	 		= [regfun;spd_reg];
		alpha 		 	= [alpha;(1e-3)*alpha];
		mref 			= [mref zeros(size(mref))];
	end
end

################################################################################################################
####### Set up inversion #######################################################################################
################################################################################################################

	### In case we use MATFREE gradMisfit to locate new RBFs
	if(locateRBFwithGrads)
		mask=0;
		Iact_free = 1.0;
		#### USE A BOUND MODEL or bounds
		a_free = 0.0;
		b_free= 1.0;

		pMis_Free = getMisfitParam(pForDip_MATFree, Wd_Dip, dobsDip, misfun, Iact_free,sback);
	end
if invertVis && invertDip
	if isSimple
		Wd   = [Wd_Dip;Wd_Vis];
		dobs = [dobsDip;dobsVis];
		pFor = [pForDip;pForVis];
		pMisRFs = getMisfitParam(pFor, Wd, dobs, misfun, Iact,sback);
	else
		modfun1 = (mAll) -> splitTwoRotationsTranslations(mAll,nDips,nShots,1);
		modfun2 = (mAll) -> splitTwoRotationsTranslations(mAll,nDips,nShots,2);
		pMisRFs1 = getMisfitParam(pForDip, Wd_Dip, dobsDip, misfun, Iact,sback,ones(length(pForDip)),modfun1);
		pMisRFs2 = getMisfitParam(pForVis, Wd_Vis, dobsVis, misfun, Iact,sback,ones(length(pForVis)),modfun2);
		pMisRFs = [pMisRFs1;pMisRFs2];
		pMisRFs1 = 0;
		pMisRFs2 = 0;
	end
	n_Moves_all = nShots + nDips;

	elseif invertDip
		pMisRFs = getMisfitParam(pForDip, Wd_Dip, dobsDip, misfun, Iact,sback);
	else
	pMisRFs = getMisfitParam(pForVis, Wd_Vis, dobsVis, misfun, Iact,sback);
end

cgit  = 50; # not used.
maxit = 5;
pcgTol = 1e-2; #not used

@everywhere function updateAngles(pFor::Array{RemoteChannel},theta_phi)
@sync begin
	@async begin
		for k=1:length(pFor)
			pFor[k] = remotecall_fetch(updateAngles,pFor[k].where,pFor[k],theta_phi);
		end
	end
end
return pFor;
end

@everywhere function updateAngles(pFor::RemoteChannel,theta_phi)
	pForMATFREE  = take!(pFor)
	worker_theta_phi = theta_phi[pForMATFREE.workerSubIdxs,:]
	pForMATFREE.theta_phi_rad=copy(worker_theta_phi);
	put!(pFor,pForMATFREE)
	return pFor;
end

function myDump(mc::Vector,Dc,iter,pInv,pMis,method="",ndips=0,noiseAngle=0,noiseTranslation=0,useVisual=false,nshots=4,plotFun=modfunForPlotting)
	# if resultsFilename!=""
		# Temp = splitext(resultsFilename);
		# Temp = string(Temp[1],"_GN",iter,Temp[2]);

			ntup = tuple((pInv.MInv.n)...);
			fullMc = plotFun(mc)[1];

		if plotting
			fullMc = reshape(IactPlot*fullMc[:],ntup);
			# fullMc = reshape(fullMc[:],ntup);
			close(888);
			figure(888);
			plotModel(fullMc,true,pInv.MInv);
			sleep(1.0);
			if(isempty(method))
				savefig(string("It",iter,"_",ntup,".png"));
				writedlm(string("Params_iter_",iter,".dat"),convert(Array{Float16},mc));
				writedlm(string("Volume_iter_",iter,"_",ntup,".dat"),convert(Array{Float16},fullMc));
			else
				visualStr ="noVis"
				if(useVisual)
					visualStr = "withVis"
				end
				savefig(string(model,"_Method_",method,"_",visualStr,"_nShots_",nshots,"_Ndips_",ndips,"_noiseAngle",noiseAngle,"_noiseTrans",noiseTranslation,"_It",iter,"_",ntup,".png"));
				writedlm(string(model,"_Method_",method,"_",visualStr,"_nShots_",nshots,"_Ndips_",ndips,"_noiseAngle",noiseAngle,"_noiseTrans",noiseTranslation,"_Params_iter_",iter,".dat"),convert(Array{Float16},mc));
				writedlm(string(model,"_Method_",method,"_",visualStr,"_nShots_",nshots,"_Ndips_",ndips,"_noiseAngle",noiseAngle,"_noiseTrans",noiseTranslation,"_Volume_iter_",iter,"_",ntup,".dat"),convert(Array{Float16},fullMc));
				if(iter > 1)
				rm(string(model,"_Method_",method,"_",visualStr,"_nShots_",nshots,"_Ndips_",ndips,"_noiseAngle",noiseAngle,"_noiseTrans",noiseTranslation,"_Volume_iter_",iter-1,"_",ntup,".dat"));
				end
			end
		end
	# end
end

pInv = getInverseParam(Mesh,modfun,regfun,alpha,mref,boundsLow,boundsHigh,
                     maxStep = 2e-1,pcgMaxIter=cgit,pcgTol=pcgTol,
					 minUpdate = 0.01, maxIter = maxit,HesPrec=HesPrec);

# end

#### Projected Gauss Newton
mc = m0;
Dc = nothing
if(method == MATBased || method == MATFree)
	pInv.maxIter = 10;
	global mc,Dc,flag,his = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
else
  pInv.maxIter = 10;
	global mc,Dc,flag,his = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
end
pInv.maxIter = 10;
myDump(mc,0,1,pInv,0,methodName,nDips,noiseAnglesDeg,noiseTrans,invertVis,nShots);

locateWithGradVis = false;

if !(method == MATBased || method == MATFree)
	outerIter = 40;
	grad_u = nothing;
	new_RBF_location = [];
	@sync for iterNum=2:outerIter
		println("Inversion using ",nRBF," basis functions")
		new_nRBF = 5;

		# Dc - the data using our method.

		u = modfunForPlotting(mc)[1];
		Ii = findall((u.>0.3) .& (u.<0.6));
		#Lets compute the gradient of the misfit function:
		if(locateRBFwithGrads)
			(m11,theta_phi1,b1) 	= splitRBFparamAndRotationsTranslations(mc,nRBF,nDips,numParamOfRBF)
			#println("sanity check: len of thetha phi:",size(theta_phi1));
			#theta_phi1 = theta_phi1[1:nDips,:];
			#println("sanity check: len of thetha phi updated:",size(theta_phi1));
			println("update angles");
			updateAngles(pForDip_MATFree,theta_phi1);
			pMis_Free = getMisfitParam(pForDip_MATFree, Wd_Dip, dobsDip, misfun, Iact_free,sback);
			println("computing grad misfit")

      if(locateWithGradVis)

      mdownsampled = zeros(size(u));
      mdownsampled=m[1:2:end,1:2:end,1:2:end]
      println("size mdown:",size(mdownsampled));
      println("size u:",size(u));
      grad_u = abs.(mdownsampled[:]-u[:]);
      grad_u = Iact_free*grad_u;
      relevantGradPoints = grad_u[Ii];

      else
			#sig = copy(u);
			grad_u = computeGradMisfit(u,Dc,pMis_Free)
			grad_u = abs.(Iact_free*grad_u[:]);
			relevantGradPoints = grad_u[Ii];
      end

		end
		if length(Ii) > new_nRBF
			if locateRBFwithGrads
				finalPoints=zeros(Int,new_nRBF);
				point=1;
				added=1;
				tries=1;
				width = 2; # determines the size of the neighborhood to ignore when adding new RBFs
				@sync while (point <=new_nRBF)
					relevantGradPoints = grad_u[Ii];
					println(point)
					added=1;

					maxVal = maximum(relevantGradPoints);
					idmax = argmax((relevantGradPoints))
					idxMaxVal = Ii[idmax];
					relevantGradPoints[idmax]=0;
					grad_u[idxMaxVal] = 0;
					#i,j,k = ind2sub(n_tup,idxMaxVal);
					(i,j,k) = Tuple(CartesianIndices(n_tup)[idxMaxVal]);
					i = max(min(i,n1-width),width+1)
					j = max(min(i,n2-width),width+1)
					k = max(min(i,n3-width),width+1)
					grad_u = reshape(grad_u,n_tup);
					grad_u[(i-width):min(i+width,n_tup[1]),(j-width):min(j+width,n_tup[2]),(k-width):min(k+width,n_tup[3])] .= 0;
					if((!(idxMaxVal in Ii ) || idxMaxVal in new_RBF_location) && (tries<length(relevantGradPoints)))
						tries = tries + 1;
						if(tries >= length(relevantGradPoints))
							grad_u = computeGradMisfit(u,Dc,pMis_Free)
							grad_u = abs.(Iact_free*grad_u[:]);
							relevantGradPoints = grad_u[Ii];
						end
						continue;
					end
					global new_RBF_location = append!(new_RBF_location,idxMaxVal);
					finalPoints[point] = idxMaxVal;
					mCenters[idxMaxVal] = 1;
					point = point + 1;
				end
				Ii = finalPoints;
				print("new RBFs added in locations:")
				println(Ii);
				print("all added RBF locations:")
				println(new_RBF_location)
			else
				Ii = Ii[randperm(length(Ii))[1:new_nRBF]];
			end
		end
		#(I1,I2,I3) = ind2sub(n_tup,I);

		(I1,I2,I3) = Tuple(CartesianIndices(n_tup)[Ii]);
		I1 = collect(Tuple(I1))[1];
		I2 = collect(Tuple(I2))[1];
		I3 = collect(Tuple(I3))[1];
		X1 = I1.*Mesh.h[1] .- 0.5*Mesh.h[1] .+ Mesh.domain[1];
		X2 = I2.*Mesh.h[2] .- 0.5*Mesh.h[2] .+ Mesh.domain[3];
		X3 = I3.*Mesh.h[3] .- 0.5*Mesh.h[3] .+ Mesh.domain[5];
		n_mc_old = numParamOfRBF*nRBF;
		n_mc_new = numParamOfRBF*nRBF + length(Ii)*numParamOfRBF;
		mc_new = zeros(length(mc) + length(Ii)*numParamOfRBF);
		if isRBF10
			mref = pInv.mref[:,1];
		end
		mref_new = zeros(Float32,length(mc) + length(Ii)*numParamOfRBF);

		mc_new[1:n_mc_old] = mc[1:n_mc_old];

		## Copying theta_phi and trans
		#(m11,theta_phi1,b1) 	= splitRBFparamAndRotationsTranslations(mc,nRBF,size(theta_phi_dip,1),numParamOfRBF)
		mc_new[(n_mc_new+1):end] = mc[(n_mc_old+1):end];
		mref_new[(n_mc_new+1):end] = mref[(n_mc_old+1):end];
		mc_new[(n_mc_old+1):numParamOfRBF:n_mc_new] .= 0.0; #ParamLevelSet.centerHeavySide - ParamLevelSet.deltaHeavySide;
		mc_new[(n_mc_old+numParamOfRBF-2):numParamOfRBF:n_mc_new] .= X1;
		mc_new[(n_mc_old+numParamOfRBF-1):numParamOfRBF:n_mc_new] .= X2;
		mc_new[(n_mc_old+numParamOfRBF):numParamOfRBF:n_mc_new] .= X3;
		pInv.boundsHigh = zeros(Float32,length(mc_new)) .+ b;
		boundsLow = zeros(Float32,length(mc_new)) .+ a;
		boundsLow[1:n_mc_old] = pInv.boundsLow[1:n_mc_old];
		if isRBF10
			mc_new[(n_mc_old+2):10:n_mc_new] .= 3.0;
			mc_new[(n_mc_old+5):10:n_mc_new] .= 3.0;
			mc_new[(n_mc_old+7):10:n_mc_new] .= 3.0;
			boundsLow[(n_mc_old+2):10:n_mc_new] .= 0.05;
			boundsLow[(n_mc_old+5):10:n_mc_new] .= 0.05;
			boundsLow[(n_mc_old+7):10:n_mc_new] .= 0.05;
		else
			mc_new[(n_mc_old+2):5:n_mc_new] .= 3.0;
			boundsLow[(n_mc_old+2):5:n_mc_new] .= 0.05;
		end
		pInv.boundsLow = boundsLow;
		mref_new[1:n_mc_new] = mc_new[1:n_mc_new];

		global mc = mc_new;
		global nRBF += new_nRBF;

		II = (sparse(1.0I, length(mc_new),length(mc_new)));
		pInv.regularizer = (m, mref, M)->TikhonovReg(m,mref,M,II);
		pInv.alpha .*= 0.8;
		if isRBF10
			spd_reg = (m,mref,M) -> RBF_SPD_regularization(m,mref,nRBF);
			pInv.regularizer  = [pInv.regularizer;spd_reg];
			pInv.mref 			  = [mref_new zeros(length(mref_new))];
		else
			pInv.mref = mref_new;
		end

		global mc, Dc,flag,his = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
		myDump(mc,0,iterNum,pInv,0,methodName,nDips,noiseAnglesDeg,noiseTrans,invertVis,nShots);

		if iterNum==outerIter
			pInv.alpha *= 0.1;
			pInv.maxIter = 10;
			mc, = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
			myDump(mc,0,iterNum+1,pInv,0,methodName,nDips,noiseAnglesDeg,noiseTrans,invertVis,nShots);
		end

	end

	if (method==RBFBasedSimple1 || method==RBF10BasedSimple1 || method==RBFBasedSimple2||method==RBF10BasedSimple2 ||method==RBF10Based ||method==RBFBased)
		meshForPlotting = getRegularMesh(Mesh.domain,Mesh.n*2);
		ntup = tuple((meshForPlotting.n)...);
		if isSimple
			modfunForPlotting = (m)->ParamLevelSetModelFunc(meshForPlotting,m;computeJacobian = 0,
				sigma = getDefaultHeavySide() ,bf = 1, numParamOfRBF = numParamOfRBF);
		else
			modfunForPlotting = (m)->ParamLevelSetModelFunc(meshForPlotting,m[1:(numParamOfRBF*nRBF)];computeJacobian = 0,
				sigma = getDefaultHeavySide() ,bf = 1, numParamOfRBF = numParamOfRBF);
		end
		fullMc = modfunForPlotting(mc)[1];
		writedlm(string(model,"_100Dips_64Shots_Noise0Deg","_Volume_final_",ntup,"_from_",tuple((Mesh.n)...),".dat"),convert(Array{Float16},fullMc));
		m = convert(Array{Float16},m);
		writedlm(string(model,"_100Dips_64Shots_Noise0Deg","_Volume_orig.dat"),m);
		writedlm(string(model,"_100Dips_64Shots_Noise0Deg","_RBF_centers.dat"),mCenters);
	end

# else
	# for k=1:40
		# sref = pInv.modelfun(mc)[1];
		# println(size(mc))
		# mref = copy(mc);
		# mref[sref.>0.5] =  1.0; #Ostu's method ?
		# mref[sref.<0.5] = 0;
		# sigma = spdiagm((sref-0.5).^2 + 1e-5);

		# tikhregfun = (m, mref, M)->TikhonovReg(m,mref,M,sigma);
		# if k==1
			# pInv.regularizer = tikhregfun#[pInv.regularizer,tikhregfun];
			# pInv.alpha       = 1e-5*(10^k)#[pInv.alpha;1e-5*(10^k)];
			# pInv.mref 		 = mref #[pInv.mref mref] #[mrefTV mref];
		# else
			 # #pInv.regularizer[2] = tikhregfun;
			 # #pInv.alpha[2]       = 1e-5*(10^k);
			 # #pInv.mref[:,2] 		= mref;
			 # pInv.regularizer = tikhregfun;
			 # pInv.alpha       = 1e-5*(10^k);
			 # pInv.mref 		= mref;
		# end
		# mc,Dc,flag,his = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit); #solveGN=projGNexplicit
		# #mc,Dc,flag,his = projGNCG(mc,pInv,pMisRFs,indCredit = [],dumpResults = dump);
		# myDump(mc,0,k+1,pInv,0,methodName,nDips,noiseAnglesDeg,noiseTrans,invertVis,nShots);
	# end
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
using NCDatasets
using Test

filename = tempname()
# The mode "c" stands for creating a new file (clobber)
ds = NCDatasets.NCDataset(filename,"c")

# define the dimension "lon" and "lat" with the size 10 and 11 resp.
NCDatasets.defDim(ds,"lon",10)
NCDatasets.defDim(ds,"lat",11)

for T in [Int32,Float32]

    local v, data

    offset = 20.
    factor = 0.1

    v = NCDatasets.defVar(ds,"scaled_var_$(T)",T,("lon","lat"), attrib = [
        "add_offset" => offset,
        "scale_factor" => factor,
    ])

    data = [-12.3*i + 23.4*j for i = 1:10, j = 1:11]

    v[:,:] = data
    @test v[:,:] ≈ data atol=1e-4

    # load without transformation (offset/scaling)
    @test v.var[:,:] ≈ (data .- offset)/factor atol=1e-4

    # write/read without transformation (offset/scaling)
    v.var[:,:] = data
    if eltype(v.var) <: Integer
        @test v.var[:,:] == round.(Int,data)
    else
        @test v.var[:,:] ≈ data
    end
end

NCDatasets.close(ds)
rm(filename)


text:
#=
This file is auto-generated. Do not edit.
=#

mutable struct InterruptibleLoad <: ControllableLoad
    name::String
    available::Bool
    bus::Bus
    model::LoadModel
    activepower::Float64
    reactivepower::Float64
    maxactivepower::Float64
    maxreactivepower::Float64
    op_cost::TwoPartCost
    _forecasts::InfrastructureSystems.Forecasts
    internal::InfrastructureSystemsInternal
end

function InterruptibleLoad(name, available, bus, model, activepower, reactivepower, maxactivepower, maxreactivepower, op_cost, _forecasts=InfrastructureSystems.Forecasts(), )
    InterruptibleLoad(name, available, bus, model, activepower, reactivepower, maxactivepower, maxreactivepower, op_cost, _forecasts, InfrastructureSystemsInternal())
end

function InterruptibleLoad(; name, available, bus, model, activepower, reactivepower, maxactivepower, maxreactivepower, op_cost, _forecasts=InfrastructureSystems.Forecasts(), )
    InterruptibleLoad(name, available, bus, model, activepower, reactivepower, maxactivepower, maxreactivepower, op_cost, _forecasts, )
end

# Constructor for demo purposes; non-functional.

function InterruptibleLoad(::Nothing)
    InterruptibleLoad(;
        name="init",
        available=false,
        bus=Bus(nothing),
        model=ConstantPower::LoadModel,
        activepower=0.0,
        reactivepower=0.0,
        maxactivepower=0.0,
        maxreactivepower=0.0,
        op_cost=TwoPartCost(nothing),
        _forecasts=InfrastructureSystems.Forecasts(),
    )
end

"""Get InterruptibleLoad name."""
get_name(value::InterruptibleLoad) = value.name
"""Get InterruptibleLoad available."""
get_available(value::InterruptibleLoad) = value.available
"""Get InterruptibleLoad bus."""
get_bus(value::InterruptibleLoad) = value.bus
"""Get InterruptibleLoad model."""
get_model(value::InterruptibleLoad) = value.model
"""Get InterruptibleLoad activepower."""
get_activepower(value::InterruptibleLoad) = value.activepower
"""Get InterruptibleLoad reactivepower."""
get_reactivepower(value::InterruptibleLoad) = value.reactivepower
"""Get InterruptibleLoad maxactivepower."""
get_maxactivepower(value::InterruptibleLoad) = value.maxactivepower
"""Get InterruptibleLoad maxreactivepower."""
get_maxreactivepower(value::InterruptibleLoad) = value.maxreactivepower
"""Get InterruptibleLoad op_cost."""
get_op_cost(value::InterruptibleLoad) = value.op_cost
"""Get InterruptibleLoad _forecasts."""
get__forecasts(value::InterruptibleLoad) = value._forecasts
"""Get InterruptibleLoad internal."""
get_internal(value::InterruptibleLoad) = value.internal