text:
"""
Data structure
--------------
mutable struct: SELF.recording

Functions
---------
SELF.recording.filter()

SELF.recording.resample()

SELF.recording.detect_spikes()
"""
module Neurons

export recording, normalize_trace, filter_trace, resample, detect_spikes, detect_spikes1, get_PSD, extract_spikes, extract_features, cluster_spikes, find_redundant, extract_isih,
		plot_detection, plot_waveforms

using DSP
using Images
using Statistics
using MultivariateStats
using GaussianMixtures
using Clustering
using LinearAlgebra
using Distances
using PyPlot

using ..MyUtils

# main type
"""
Fields
------
+ file::String
+ sf::Real
+ sf_preprocessed::Real
+ n_channels::Real
+ trace::AbstractArray # supposed to contain sample X channel
+ trial_marker::AbstractArray
+ trace_preprocessed::AbstractArray
+ preprocessing_steps::Vector{String}
+ spike_idx::AbstractArray
+ detection_threshold::AbstractArray
+ isih::AbstractArray
+ spikeforms::AbstractArray # channel (val_source) X channel (t_source) X time of spike/spikeform X sample
+ features::AbstractArray # channel X feature space representation of spikes
+ spike_class::AbstractArray # channel X clusters (as spike indices)
+ BIC::AbstractArray
"""
mutable struct recording
    ### FIELDS
    file::String
    sf::Real
    sf_preprocessed::Real
    n_channels::Real
    trace::Array{Array{Float64,2}}#AbstractArray # supposed to contain sample X channel
    trial_marker::AbstractArray
    trace_preprocessed::Array{Array{Float64,2}}#AbstractArray
    preprocessing_steps::Vector{String}
    spike_idx::AbstractArray
	detection_threshold::AbstractArray
    isih::AbstractArray
    spikeforms::AbstractArray # channel (val_source) X channel (t_source) X time of spike/spikeform X sample
    features::AbstractArray # channel X feature space representation of spikes
    spike_class::AbstractArray # channel X clusters (as spike indices)
    BIC::AbstractArray

    ### CONSTRUCTOR
    # recording(file, sf, trace, trial_marker) =
    #    new(file, sf, 0, trace, trial_marker, Array{Float64,2}(undef,0,0), Vector{String}(undef,0), Array{Float64,3}Array{Float64,2}(undef,0,0,0))
    # recording() = new("", 0, 0, Array{Float64,2}(undef,0,0), Array{Float64,2}(undef,0,0), Array{Float64,2}(undef,0,0), Vector{String}(undef,0), Array{Float64,3}(undef,0,0,0))
    recording() = new("", 0, 0, 0, [], [], [], Vector{String}(undef,0), [], [], [], [], [], [])
end

# functions

function normalize_trace(n::Neurons.recording, what="trace")

	x = deepcopy( getfield( n, Symbol(what) ) )
	sd = std.(x)
	mn = mean.(x)
	[x[ii]=(x[ii] .- mn[ii]) ./ sd[ii] for ii in 1:n.n_channels] ### normalization
	n.trace_preprocessed = x
	n.preprocessing_steps = ["normalized"];

end

"""
Documentation ...
"""
function filter_trace(n::Neurons.recording, filter_type::String, freqs::AbstractArray, filter_design::AbstractArray, filter_order::AbstractArray, bandwidth::AbstractArray, from_chnnls::AbstractArray, what::String="trace")
"""
    filter(n::recording, band::Tuple{Real}, filter_type)
    filter_type: lowpass, highpass, bandpass, notch
    freqs = highpass: freq | lowpass: freq | passband (Wn1,Wn2) frequency | Notch: array of frequencies | CommonAverage: channels used for averaging
    what: spcifies whether filter works on trace or trace_preprocessed. If the latter trace_preprocessed is replaced, respectively.
"""

    ######################
    ### raw trace
    ######################
    if what == "trace"
        nyq = n.sf/2;
        freqs_n = freqs./nyq; # normalized to nquist frequncy
        n.trace_preprocessed = []
        n.sf_preprocessed = n.sf # restore original sampling frequency

        if filter_type == "Bandpass"
            rt = Bandpass( freqs_n[1], freqs_n[2] );
            dm = getfield( Main, Symbol(filter_design[1]) )(filter_order[1]);
            f = digitalfilter(rt, dm);
            ### loop over channels and filter trace
            for trc in n.trace
                #trc_padded = flatten([[.0 for ii in 1:pads], trc, [.0 for ii in 1:pads]])
                #push!(n.trace_preprocessed, filtfilt(f, trc_padded)[pads+1:end-pads]);
                push!(n.trace_preprocessed, filtfilt(f, trc));
            end
        elseif filter_type == "Lowpass"
            rt = Lowpass( freqs_n[1] );
            dm = getfield( Main,Symbol(filter_design[1]) )(filter_order[1]);
            f = digitalfilter(rt, dm);
            for trc in n.trace
                # trc_padded = flatten([[.0 for ii in 1:pads], trc, [.0 for ii in 1:pads]])
                # push!(n.trace_preprocessed, filtfilt(f, trc_padded)[pads+1:end-pads]);
                push!(n.trace_preprocessed, filtfilt(f, trc));
            end
        elseif filter_type == "Highpass"
            rt = Highpass( freqs_n[1] );
            dm = getfield( Main,Symbol(filter_design[1]) )(filter_order[1]);
            f = digitalfilter(rt, dm);
            for trc in n.trace
                # trc_padded = flatten([[.0 for ii in 1:pads], trc, [.0 for ii in 1:pads]])
                # push!(n.trace_preprocessed, filt(f, trc_padded)[pads+1:end-pads]);
                push!(n.trace_preprocessed, filtfilt(f, trc));
            end
        elseif filter_type == "Notch"
            for trc in n.trace
                trc = flatten(trc)
                for ff in freqs_n.*nyq
                    nn = iirnotch(ff, bandwidth[1]; fs=n.sf)
                    trc = filtfilt(nn, trc);
                end
                push!(n.trace_preprocessed, trc);
            end
        elseif filter_type == "CommonAverage"
            chnnls = freqs
            average_trace = [median([n.trace[ch][idx] for ch in chnnls]) for idx in 1:length(n.trace[1])]
            # n.trace_preprocessed = [channel .- average_trace for channel in n.trace[from_chnnls]]
			[n.trace_preprocessed[ch] = n.trace_preprocessed[ch] .- average_trace for ch in from_chnnls]
            #n.trace_preprocessed = [ch .- mean(n.trace) for ch in n.trace]
        end
        n.preprocessing_steps = [string( filter_type, ": ",freqs )];

    ######################
    ### preprocessed trace
    ######################
    elseif what == "trace_preprocessed" && !isempty(n.preprocessing_steps)
        nyq = n.sf_preprocessed/2;
        freqs_n = freqs./nyq; # normalized to nquist frequncy

        if filter_type == "Bandpass"
            rt = Bandpass( freqs_n[1], freqs_n[2] );
            dm = getfield( Main,Symbol(filter_design[1]) )(filter_order[1]);
            f = digitalfilter(rt, dm);
            for (itrc, trc) in enumerate(n.trace_preprocessed)
                # trc_padded = flatten([[.0 for ii in 1:pads], trc, [.0 for ii in 1:pads]])
                # n.trace_preprocessed[itrc] = filt(f, trc_padded)[pads+1:end-pads];
                n.trace_preprocessed[itrc] = filtfilt(f, trc);
            end
        elseif  filter_type == "Lowpass"
            rt = Lowpass( freqs_n[1] );
            dm = getfield( Main,Symbol(filter_design[1]) )(filter_order[1]);
            f = digitalfilter(rt, dm);
            for (itrc, trc) in enumerate(n.trace_preprocessed)
                # trc_padded = flatten([[.0 for ii in 1:pads], trc, [.0 for ii in 1:pads]])
                # n.trace_preprocessed[itrc] = filt(f, trc_padded)[pads+1:end-pads];
                n.trace_preprocessed[itrc] = filtfilt(f, trc);
            end
        elseif filter_type == "Highpass"
            rt = Highpass( freqs_n[1] );
            dm = getfield( Main,Symbol(filter_design[1]) )(filter_order[1]);
            f = digitalfilter(rt, dm);
            for (itrc, trc) in enumerate(n.trace_preprocessed)
                # trc_padded = flatten([[.0 for ii in 1:pads], trc, [.0 for ii in 1:pads]])
                # n.trace_preprocessed[itrc] = filt(f, trc_padded)[pads+1:end-pads];
                n.trace_preprocessed[itrc] = filtfilt(f, trc);
            end
        elseif filter_type == "Notch"
            for (itrc, trc) in enumerate(n.trace_preprocessed)
                trc = flatten(trc)
                for ff in freqs_n.*nyq
                    nn = iirnotch(ff, bandwidth[1]; fs=n.sf_preprocessed)
                    trc = filtfilt(nn, trc);
                end
                n.trace_preprocessed[itrc] = trc;
            end
        elseif filter_type == "CommonAverage"
            chnnls = freqs
            average_trace = [median([n.trace_preprocessed[ch][idx] for ch in chnnls]) for idx in 1:length(n.trace_preprocessed[1])]
            # n.trace_preprocessed = [channel .- average_trace for channel in n.trace_preprocessed]
			[n.trace_preprocessed[ch] = n.trace_preprocessed[ch] .- average_trace for ch in from_chnnls]
			#n.trace_preprocessed = [ch .- mean(n.trace_preprocessed) for ch in n.trace_preprocessed]
        end
        push!(n.preprocessing_steps, string( filter_type, ": ", freqs ));
    else
        print(n.preprocessing_steps)
        @warn "SELF.trace_preprocessed is empty. Process SELF.trace first"
    end
end

"""
Arguments:

+ n::Neurons.recording
+ chans::Array{Int64,1}
+ SNR::Array{Float64,1}
+ lockout::Int64=30;
+ polarity::String="negative",
+ ww::Array{Int64,1}=[5,15],
+ stp::Int64=1,
+ w_lockout::Array{Int64,1}=[10],
+ method::String="petrantonakis"

"""
function detect_spikes(
	n::Neurons.recording,
	chans::Array{Int64,1},
	SNR::Array{Float64,1},
	lockout::Int64=30;
	#
	polarity::String="negative",
	ww::Array{Int64,1}=[5,15],
	stp::Int64=1,
	w_lockout::Array{Int64,1}=[10],
	method::String="petrantonakis")
#
#   spikes are stored as:s
#   n.spike_idx: channel X w X Dw/indices of peaks (Dw only if method is "petrantonakis")
#
#   w: window length if method is "petrantonakis"
#   (for consistency if method is "threshold" the same structure is stored )

### WAVEFORM/THRESHOLD BASED SPIKE DETECITON
    if method == "petrantonakis" ### ref: https://www.frontiersin.org/articles/10.3389/fnins.2015.00452/full
        # Dw = [ [[] for ii in 1:n.n_channels] for jj in eachindex(ww) ] # ww X channel
        # Spks = [ [CartesianIndex{1}[] for ii in 1:n.n_channels] for jj in eachindex(ww) ] # ww X channel
        Dw_spk = [ [[[],[]] for ii in eachindex(ww)] for jj in 1:n.n_channels ] ### ch X w X Dw,spk_idx
		n.detection_threshold = []

        # for (ii,ch) in enumerate([n.trace_preprocessed[2]]) ### loop over channels
		for ii in chans ### loop over channels

			ch = n.trace_preprocessed[ii]

			for (jj,w) in enumerate(ww) ### loop over w

                ### compute the Dw signal (euklidean distance signal)
                f = ch[collect(0:w-1).+collect(1:stp:length(ch)-w+1)'] ### cut signal in overlapping chunks
				Dw_spk[ii][jj][1] = flatten(sqrt.(sum(diff(f, dims=2).^2, dims=1))); ### compute euclidean distance between consecutive chunks

                ### normalization of Dw
				sig = median( Dw_spk[ii][jj][1] ./ 0.6745 )
				thr = SNR[jj]*sig
				append!(n.detection_threshold, thr)
                # Dw_spk[ii][jj][1] = Dw_spk[ii][jj][1] ./ th ### store the standardized Dw signal
                ### find all local mixima in Dw
                loc_maxima = findlocalmaxima(Float64.(Dw_spk[ii][jj][1])) ### find all local maxima
                ### remove local mixima < thr
                Dw_spk[ii][jj][2] = loc_maxima[findall(Dw_spk[ii][jj][1][loc_maxima] .> thr)] ### extract all local maxima > SNR, SNR depends on w

                ### remove redundant triggered spikes and keep the maximum of each
                if length(Dw_spk[ii][jj][2]) > 1

					Dw_spk[ii][jj][2] = find_redundant( Dw_spk[ii][jj][1], Dw_spk[ii][jj][2], lockout );

				end

			end #for w: jj

            ### remove spikes that are detected only by one value of w
			loop_cntr = 0
			while length(Dw_spk[ii][1][2]) != length(Dw_spk[ii][2][2]) ### necessary if there a spikes in value of w with less spikes that are not in w with more spikes

				loop_cntr += 1

				if loop_cntr == 10 ### 2 should be enough

					break;

				end

				if !any( [(length(Dw_spk[ii][1][2])==0), (length(Dw_spk[ii][2][2])==0)] )

					if length(Dw_spk[ii][1][2]) > length(Dw_spk[ii][2][2])

						dd = pairwise(Euclidean(), cartesian2matrix(Dw_spk[ii][2][2])', cartesian2matrix(Dw_spk[ii][1][2])', dims=2)
					    deleteat!(Dw_spk[ii][1][2], findall(flatten(.!any(dd .< w_lockout, dims=1)))) ### delete all spikes in Dw_spk[ii][1][2] that do not have a partner spike in Dw_spk[ii][2][2]

					elseif length(Dw_spk[ii][1][2]) < length(Dw_spk[ii][2][2])

						dd = pairwise(Euclidean(), cartesian2matrix(Dw_spk[ii][1][2])', cartesian2matrix(Dw_spk[ii][2][2])', dims=2)
					    deleteat!(Dw_spk[ii][2][2], findall(flatten(.!any(dd .< w_lockout, dims=1))))

					end

				else ### if one has no spikes at all

					Dw_spk[ii][1][2] = [];
					Dw_spk[ii][2][2] = [];
					# Dw_spk[ii][1][1] = []
					# Dw_spk[ii][2][1] = []

				end

			end

			if length(Dw_spk[ii][1][2]) != length(Dw_spk[ii][2][2])

				@info "@Channel ", ii, ": different spike numbers for values of w"

			end

			### scatter plot of w1 X w2
	        if (length(spk_idx_w1)-length(spk_idx_w1) == 0) && length(recording.spike_idx[channel][1][2]) != 0

	            ax1 = subplot(1,3,1)
	            w1 = recording.spike_idx[channel][1][1][spk_idx_w1]
	            w2 = recording.spike_idx[channel][2][1][spk_idx_w2]
	            pols = cartesian2polar(w1, w2)
	            r = pols[:,1]
	            phi = pols[:,2]
	            ax1.scatter( r, phi, color=pcols[2])

	        end

        end #for channels: ii

    ### THRESHOLD BASED SPIKE DETECTION
    elseif method == "threshold"

		Dw_spk = [ [[[],[]] for ii in [1]] for jj in 1:n.n_channels ] ### ch X w X DW,spk_idx
		n.detection_threshold = ones(n.n_channels) .* NaN
        x = n.trace_preprocessed[chans]

        sig = [median(abs.(ch) ./ 0.6745 ) for ch in x]
		thr = SNR[1] .* sig
		n.detection_threshold[chans] = thr

        # x = .-x ./ th

        if polarity == "both"

			ind_max = findlocalmaxima.(x)
            ind_max = [[lm[1] for lm in ch] for ch in ind_max] # cartesian coordinates to indices
            ind_min = findlocalminima.(x)
            ind_min = [[lm[1] for lm in ch] for ch in ind_min]
            ind_max = [ind_max[ii][findall(x[ii][ind_max[ii]] .> thr[ii])] for ii in eachindex(chans)]    # local minima with minimal SNR
            ind_min = [ind_min[ii][findall(x[ii][ind_min[ii]] .< -thr[ii])] for ii in eachindex(chans)]    # local maxima with minimal SNR
            ind = [sort(flatten([ind_max[ii], ind_min[ii]])) for ii in 1:length(ind_max)] # channel X index_of_spikes

		elseif polarity == "positive"

			ind = findlocalmaxima.(x)
            ind = [[lm[1] for lm in ch] for ch in ind] # cartesian coordinates to indices
            ind = [ind[ii][findall(x[ii][ind[ii]] .> thr[ii])] for ii in eachindex(chans)]    # local minima with minimal SNR

		elseif polarity == "negative"

			ind = findlocalminima.(x)
            ind = [[lm[1] for lm in ch] for ch in ind] # cartesian coordinates to indices
            ind = [ind[ii][findall(x[ii][ind[ii]] .< -thr[ii])] for ii in eachindex(chans)]

		end

        ### remove redundantly triggered spikes and store the maximum respectively
        for (ii,ch) in enumerate(chans)

			chans_spike_idxs = ind[ii]
    		Dw_spk[ch][1][2] = CartesianIndex.( find_redundant( x[ii], chans_spike_idxs, lockout ) )

		end

    end # method

    ### store clean indices in struct
    # n.spike_idx = ind_clean#[ind_clean[ii][d[ii]] for ii in 1:length(ind_clean)]
    n.spike_idx = Dw_spk; # channel X w X Dw/index of peak (Dw only if method is "petrantonakis")
end

"""
Documentation ...
"""
function extract_isih(n::Neurons.recording)
### Extracts the inter-spike intervals
    dt = 1/n.sf_preprocessed
    isih = []
    for ch in  eachindex(n.spike_idx)
        push!(isih, [])
        for w in eachindex(n.spike_idx[1])
            push!(isih[ch], [])
            if length(n.spike_idx[ch][w][2]) > 1
                isih[ch][w] = cartesian2matrix(diff(n.spike_idx[ch][w][2])) .* dt
            end
        end
    end
    n.isih = isih # value_source X w X time_source X time/values
end

"""
Documentation ...
"""
function extract_spikes(n::Neurons.recording, mp, w::Int64=1)
### Extract spike spikeforms and time stamps in ms.
#   extracts the spikeforms at indices n.spike_idx (normalized around zero)
#   from the filtered signal x using a fixed window around the
#   window +-mp of the spikes.
#	ONLY FOR ONE VALUE OF w
#   structure: n.spikeforms[val_source][time_source][spike_time/value]
#

    dt = 1/n.sf_preprocessed
	window = collect(mp[1]:1:mp[2]);
	spkfrms = []

	for ii in 1:n.n_channels # loop over channels as value source

		push!(spkfrms, [])
	    # for w in 1:size(n.spike_idx[1],1)
	    #     push!(spkfrms[ii], [])

			for (jj,spk_idxs) in enumerate(n.spike_idx) # loop over channels for time source

				push!(spkfrms[ii], [[],[]]) ### t, s

				for spk in spk_idxs[w][2] # loop over individual spikes

				    # skip spikes that have negative indices
	                spike_window_idxs = spk.I .+ window # generate window for each spike
	                spike_window_times = spike_window_idxs .* dt .* 1000 # spike times in ms

					if spike_window_idxs[1] > 0 && spike_window_idxs[end] <= length(n.trace_preprocessed[ii])# skip spikes that are not complete in window

						push!( spkfrms[ii][jj][1], spike_window_times ) # store spike times
	                    push!( spkfrms[ii][jj][2], n.trace_preprocessed[ii][spike_window_idxs] .- mean( n.trace_preprocessed[ii][spike_window_idxs]) ) #store spike values

					end
	            end
	        end
	    # end
	end
	n.spikeforms = spkfrms # value_source X w X time_source X time/values

end

"""
Documentation ...
"""
function extract_features(n::Neurons.recording)
    n.features = []
    for ch in 1:n.n_channels
        push!(n.features, [])
        if !isempty(n.spikeforms[ch][ch][2])
            ### make matrix for PCA
            m = zeros(length(n.spikeforms[ch][ch][2]), length(n.spikeforms[ch][ch][2][1])); # #spikes X #samples per spike
            [m[ii,:] = n.spikeforms[ch][ch][2][ii] for ii = 1:length(n.spikeforms[ch][ch][2])]; # spike X sample
            ### fit PCA
            M = fit(PCA, m)
            n.features[ch] = projection(M) # projection(M): spike X PC
        end
    end
end

"""
Documentation ...
"""
function cluster_spikes(n::Neurons.recording, algorithm::AbstractArray, init_method::AbstractArray, feats::AbstractArray, n_clusters::AbstractArray, nIter_EM::AbstractArray, nIter_kmeans_init::AbstractArray, cov_kind::AbstractArray)

############
# cluster given features using _algorithm_
# BIC optimization for cluster number
#
########
# algorithm: kmeans
#
########
# algorithm: gmm
# init_method: :split or :kmeans

    n.spike_class = []
    ###
    if algorithm[1] == "kmeans"
        for ch in 1:n.n_channels
            #push!(n.spike_class, [])
            if size(n.features[ch], 2) > feats[end] # if there are enough features
                result = kmeans(n.features[ch][:,feats]', n_clusters[1]);
                clusters = []
                for cl = 1:n_clusters[1]
                    push!(clusters, findall(result.assignments .== cl))
                end
                push!(n.spike_class, clusters)
            else ### push empty array for that channel
                push!(n.spike_class, [])
            end
        end
    ###
    elseif algorithm[1] == "gmm"
        n.BIC = []
        ### for each channel fit gmm with n components specified in
        ### n_clusters, compute BIC fear each value, find minimum BIC and use
        ### for clustering
        for ch in 1:n.n_channels
            if size(n.features[ch],2) > length(feats) ### if enough features
                N = size(n.features[ch][:,feats], 1)
                BICS = Array{Float64,2}(undef, 0, 2)
                for cmpnts in n_clusters
                    if size(n.features[ch], 1) > cmpnts
                        M = GMM(cmpnts, n.features[ch][:,feats];
                            method=init_method[1], kind=cov_kind[1], nInit=nIter_kmeans_init[1], nIter=nIter_EM[1], nFinal=nIter_EM[1])
                        LL = sum(llpg(M, n.features[ch][:,feats]))
                        par = cmpnts + cmpnts*length(feats) + cmpnts* length(feats)*(length(feats)+1)/2 # number of parameters
                        bic = -2*LL + log.(N)*par
                        BICS = vcat(BICS, [cmpnts bic])
                    end
                end
                push!(n.BIC, BICS)
                ### fit model with optimal number of components and cluster spikes
                opt_cmpnts = Int(BICS[findmin(BICS[:,2])[2]])
                push!(n.spike_class, [])
                M = GMM(opt_cmpnts, n.features[ch][:,feats];
                    method=init_method[1], kind=cov_kind[1], nInit=nIter_kmeans_init[1], nIter=nIter_EM[1], nFinal=nIter_EM[1])
                posteriors = gmmposterior( M, n.features[ch][:,feats] )
                classes = [findmax(posteriors[1], dims=2)[2][spk][2] for spk in 1:length(findmax(posteriors[1], dims=2)[2])] ### extract cluster index for each spike from CartesianCoordinates
                n.spike_class[ch] = [[] for ii in 1:opt_cmpnts]
                for (spk_idx,cl) in enumerate(classes)
                    push!(n.spike_class[ch][cl], spk_idx)
                end
            else
                push!(n.spike_class, [])
                push!(n.BIC, [])
            end
        end
        # for ch in 1:n.n_channels
        #     print(ch)
        #     push!(n.spike_class, [])
        #     if size(n.features[ch],2) > length(feats) ### if enough features
        #         M = GMM(n_clusters[1], n.features[ch][:,feats];
        #             method=:kmeans, kind=cov_kind[1], nInit=nIter_kmeans_init[1], nIter=nIter_EM[1], nFinal=nIter_EM[1])
        #         posteriors = gmmposterior( M, n.features[ch][:,feats] )
        #         classes = [findmax(posteriors[1], dims=2)[2][spk][2] for spk in 1:length(findmax(posteriors[1], dims=2)[2])] ### extract cluster index for each spike from CartesianCoordinates
        #         n.spike_class[ch] = [[] for ii in 1:n_clusters[1]]
        #         for (spk_idx,cl) in enumerate(classes)
        #             push!(n.spike_class[ch][cl], spk_idx)
        #         end
        #     else
        #         n.spike_class[ch] = [[] for ii in 1:n_clusters[1]]
        #     end
        # end
    end
end

"""
Up- or downsampling
Acts on either n.trace or n.trace_preprocessed
"""
function resample(n::Neurons.recording, frequency::Real)
    #idx = 1:factor:size(n.trace)[1]
    #n.trace_preprocessed = n.trace[idx];
    #n.sf_decimated = n.sf/factor;
    if isempty(n.preprocessing_steps)
        n.trace_preprocessed = DSP.resample( n.trace,frequency//convert(Int64, n.sf) );
    else
        n.trace_preprocessed = DSP.resample( n.trace_preprocessed,frequency//convert(Int64, n.sf_preprocessed) );
    end
    n.sf_preprocessed = frequency;
    push!(n.preprocessing_steps, string( "Decimated to: ",frequency ))
end

"""
Up- or downsampling
Acts on either n.trace or n.trace_preprocessed
"""
function get_PSD(n::Neurons.recording, which, psd_type)
    if which == "trace"
        return [getfield( Main, Symbol(psd_type) )( flatten(ch), fs=n.sf ) for ch in n.trace]
    elseif which == "trace_preprocessed"
        return [getfield( Main, Symbol(psd_type) )( flatten(ch), fs=n.sf_preprocessed ) for ch in n.trace_preprocessed]
    end
end

function find_redundant(ref::AbstractArray, ind::AbstractArray, lockout::Int)
    ### ref: 1D array containing amplitude like values
    ### ind: 1D array containing indices of identified events, e.g. spikes
    ### lockout: Integer defining minimum distance of evvents in samples that are treated as separate events

    ### make sure it's cartesian!
    ind = CartesianIndex.(ind)
    ### matrix with all pairs of spike indices
    # mm = Iterators.product(flatten(cartesian2matrix(ind)), flatten(cartesian2matrix(ind)))
	### differences between all pairs of spike indices
	# dd = [ float(abs(diff([x for x in m])[1][1])) for m in mm ]
	dd = pairwise(Euclidean(), flatten(cartesian2matrix(ind))', flatten(cartesian2matrix(ind))', dims=2)
	dd[diagind(dd)] .= NaN # set diagonal to NaN
	sup_diag = dd[diagind(dd, 1)]
	gg = findall(x -> x > lockout, sup_diag )

	if all(sup_diag .> lockout) ### no double triggered spikes

		groups = [[spk_idx] for spk_idx in ind]

	elseif all(sup_diag .< lockout) ### only one double triggered spike

		groups = [ind]

	else ### single and double triggered spikes

		groups = []
	    current_group = []

		for (ij,spk_id) in enumerate(ind)

			push!(current_group, spk_id)

			if in(ij, gg) || ij == length(ind) ### ij == length(ind) required for last element

				push!(groups, current_group)
	            current_group = []

			end
		end
	end
	### find maxima of groups and store indices
    ind_clean = []
    for gr in groups

		max_id = findmax( abs.(ref[CartesianIndex.( gr )]) )
		# max_id = findmax(ref[gr])
    	push!(ind_clean, gr[max_id[2]])

	end

    #################
    ### DELETE ME ###a
    #################
    # ind_clean = []
    # ind_ = deepcopy(ind)
    # if length(ind_) > 1
    #     groups = []
    #     current_group = []
    #     push!(current_group, ind_[1])
    #     deleteat!(ind_, 1)
    #
    #     ### extract groups of detected event that are too close to each other to represent seperate spikes
    #     while !isempty(ind_)
    #         ### while next data point is to close to preceeding one, add it to current group
    #         while ind_[1][1]-current_group[end][1] < lockout ### double index because CartesianIndex
    #             push!(current_group, ind_[1])
    #             deleteat!(ind_, 1)
    #             if isempty(ind_) ### if already at the last spike index
    #                 break;
    #             end
    #         end
    #
    #         ### add current_group to group
    #         push!(groups, current_group)
    #         current_group = []
    #         if !isempty(ind_) ### deals with if last index was already deleted in while loop
    #             push!(current_group, ind_[1])
    #             deleteat!(ind_, 1)
    #         end
    #     end
    #
    #     ### for each group extract the maximum as representation
    #     ind_c = []
    #     for gr in groups
    #         if length(gr) > 1 ### if the group is a "group"
    #             gg = [gr[ii] for ii in eachindex(gr)]
    #             push!( ind_c, gg[findmax(ref[gg])[2]] )
    #         else ### if the group has only one entry
    #             push!(ind_c, gr[1])
    #         end
    #     end
    #     push!(ind_clean, ind_c)
    # else ### if no spikes in channel, add empty list to kepp size consistent
    #     push!(ind_clean, [])
    # end

    return ind_clean
end

################
### PLOTTING ###
################

function plot_detection(
        recording::Neurons.recording;
        channel::Int64,
        method::String="petrantonakis",
        polarity::String="negative")

	colors = myColors() # from MyHelpers.jl
	pcols = map(col -> (red(col), green(col), blue(col)), colors)

    if method == "petrantonakis"

        fig = figure(figsize=(13,5))
        rc("xtick",labelsize=8)
        rc("ytick",labelsize=8)

		### time course of Dw signal
        Dw_w1 = recording.spike_idx[channel][1][1]
        Dw_w2 = recording.spike_idx[channel][2][1]
		### indices of spikes
		spk_idx_w1 = recording.spike_idx[channel][1][2]
        spk_idx_w2 = recording.spike_idx[channel][2][2]
		### spike indices from cartesian to matrix
		if !isempty(spk_idx_w1) && !isempty(spk_idx_w2)
			spk_idx_w1 = cartesian2matrix(spk_idx_w1)
			spk_idx_w2 = cartesian2matrix(spk_idx_w2)
		end
        w1_threshold = recording.detection_threshold[1]
        w2_threshold = recording.detection_threshold[2]

        ### scatter plot of w1 X w2
        if (length(spk_idx_w1)-length(spk_idx_w1) == 0) && length(recording.spike_idx[channel][1][2]) != 0

            ax1 = subplot(1,3,1)
            w1 = recording.spike_idx[channel][1][1][spk_idx_w1]
            w2 = recording.spike_idx[channel][2][1][spk_idx_w2]
            pols = cartesian2polar(w1, w2)
            r = pols[:,1]
            phi = pols[:,2]
            ax1.scatter( r, phi, color=pcols[2])

        end

        ax2 = subplot(1,3,2)
        ax2.plot( Dw_w1, color=pcols[1], zorder=1 )
        ax2.plot( 1:length(Dw_w1), ones(length(Dw_w1))*w1_threshold, color=:red, linestyle="--" )
        if length(recording.spike_idx[channel][1][2]) > 0

            spk_idx_w1 = cartesian2matrix(recording.spike_idx[channel][1][2])
            ax2.scatter(spk_idx_w1, Dw_w1[spk_idx_w1], color=pcols[2], zorder=2)

        end

        ax3 = subplot(1,3,3)
        ax3.plot(Dw_w2, color=pcols[1], zorder=1)
        ax3.plot( 1:length(Dw_w2), ones(length(Dw_w2))*w2_threshold, color=:red, linestyle="--" )
		if length(recording.spike_idx[channel][2][2]) > 0

            spk_idx_w2 = cartesian2matrix(recording.spike_idx[channel][2][2])
            ax3.scatter(spk_idx_w2, Dw_w2[spk_idx_w2], color=pcols[2], zorder=2)

        end

    elseif method == "threshold"

        fig = figure(figsize=(13,5))
        rc("xtick",labelsize=8)
        rc("ytick",labelsize=8)

        x = recording.trace_preprocessed[channel]
        spk_idx = cartesian2matrix(recording.spike_idx[channel][1][2])

		if polarity == "negative"

			thr = -recording.detection_threshold[channel]

		elseif polarity == "positive"

			thr = recording.detection_threshold[channel]

		elseif polarity == "both"

			thr = [-recording.detection_threshold[channel], recording.detection_threshold[channel]]
        end

        plot( x )
        scatter( spk_idx, x[spk_idx], color=pcols[2], zorder=2 )

        if length(thr) != 2

			plot( 1:length(x), ones(length(x))*thr, color=:red, linestyle="--" )

		else

			plot( 1:length(x), ones(length(x))*thr[1], color=:red, linestyle="--" )
            plot( 1:length(x), ones(length(x))*thr[2], color=:red, linestyle="--" )
        end
    end
end

function plot_waveforms(
	recording::Neurons.recording,
	t::Array{Float64,1},
	ylims::Array{Int64,2},
	w::Int64=1;
	#
	t_source )
	##########################
	##########################
	##########################

	colors = myColors() # from MyHelpers.jl
	pcols = map(col -> (red(col), green(col), blue(col)), colors)
	half_window = abs(minimum(t))/1000

	fig = figure(figsize=(15,7), dpi=300)
	rc("xtick",labelsize=4)
	rc("ytick",labelsize=4)

	# ylims = [-1000 1000]

	if t_source == "val_source"
		corresponding = true
	else
		corresponding = false
	end

	for val_source = 1:recording.n_channels
	    #### t_cource determines whose channels spike times are used to cut out the the waveforms: if t_cource = val_source then channel and spiketimes are matched, i.e. 1<->1, 2<->2 etc
		if corresponding
			t_source = val_source
		end
	#     t_source = val_source

	    ### sort spikes according to amplitude
	    sorted_idx = sortperm([maximum(spk) for spk in recording.spikeforms[val_source][t_source][2]] .- [minimum(spk) for spk in recording.spikeforms[val_source][t_source][2]], rev=true )
	    ### how many spikes to plot
	    # n = length(sorted_idx)
	    n = 10

	    if !isempty(sorted_idx) ### if there are spikes in sorted_idx, i.e. if channel has spikes

	        spks_to_plot = []

	        ### loop over clusters of current channel
	        if !isempty(recording.spike_class)
	            for cl in recording.spike_class[t_source] ### plot the clusters based on the t_source channel
	                push!(spks_to_plot, intersect(cl, sorted_idx))
	            end
	        else
	            push!(spks_to_plot, sorted_idx) ### i.e. one "cluster" if not yet clustered
	        end

	        ### check if there are enough spikes if number of spikes to plot is specified
	        if length(sorted_idx) < n
	            n = length(sorted_idx)
	        end

	        ### plot spikes belonging to clusters
			offset = 0
	        for (ii,spk_cl) in enumerate(spks_to_plot)
	            if !isempty(spk_cl) ### cluster is not empty
	                subplot(2,5,val_source)
				   ### plot individual traces
				   plot(t, reduce(hcat, recording.spikeforms[val_source][t_source][2][spk_cl]),
					  color=pcols[ii], alpha=0.3, linewidth=0.1, zorder=1) ### plot individual spikes
				   ### plot mean waveform +- std
				   clusters_mean = mean(reduce(hcat, recording.spikeforms[val_source][t_source][2][spk_cl]), dims=2) .+ offset
				   clusters_std = std(reduce(hcat, recording.spikeforms[val_source][t_source][2][spk_cl]), dims=2)
				   plot(t, clusters_mean, color=pcols[ii], alpha=0.5, linewidth=0.5, zorder=2) ### plot mean spike waveform of cluster
				   fill_between(t, flatten(clusters_mean+clusters_std), flatten(clusters_mean-clusters_std), color=pcols[ii], alpha=0.2, linewidth= 0.1)
				   ###
				   title(string("Channel ", val_source, " (", n, ")"), fontsize=5)
				   ylim(ylims[1], ylims[2])
	# 			   legend()
	            end
				offset = offset +0
	        end

	        ### add x- and ylabel
	        if val_source in collect(6:10)
	            subplot(2,5,val_source)
	            xlabel("ms", fontsize=5)
	        elseif val_source in [1,6]
	            subplot(2,5,val_source)
	#             ylabel("amplitude")
	        end
			if val_source in collect(1:5)
			   xticks(-half_window*1000:1:half_window*1000, []) ### half_window comes from waveform extraction step
			end
			if val_source in [2:5 7:10]
			   yticks(ylims[1]:100:ylims[2], [])
			end

	    else ### plot empty axes of channel has no spikes
	        subplot(2,5,val_source)
	        plot(t, t,
	            linewidth=0)
	        title(string("Channel ", val_source, " (", n, ")"))
	        ### add x- and ylabel
	        ### add x- and ylabel
	        if val_source in collect(6:10)
	            subplot(2,5,val_source)
	            xlabel("ms", fontsize=5)
	        elseif val_source in [1,6]
	            subplot(2,5,val_source)
	#             ylabel("amplitude")
	        end
			if val_source in collect(1:5)
			   xticks(-4:2:4, [])
			end
			if val_source in [2:5 7:10]
			   yticks(ylims[1]:100:ylims[2], [])
			end

	    end

	end

end

###
end

max_stars_repo_name: maalaria/Newbe

max_issues_repo_name: maalaria/Newbe

max_forks_repo_name: maalaria/Newbe


text:
using RobotDynamics
using LinearAlgebra
using StaticArrays
using Rotations
using BenchmarkTools
using Test

struct SatellitePair{R,T} <: LieGroupModel
    J1::SMatrix{3,3,T,9}   # inertia of satellite 1
    J2::SMatrix{3,3,T,9}   # inertia of satellite 2
end

function RobotDynamics.dynamics(model::SatellitePair, x, u)
    vs = RobotDynamics.vec_states(model, x)
    qs = RobotDynamics.rot_states(model, x)
    ω1 = vs[2]
    ω2 = vs[3]
    q1 = qs[1]
    q2 = qs[2]

    J1, J2 = model.J1, model.J2
    u1 = u[SA[1,2,3]]
    u2 = u[SA[4,5,6]]
    ω1dot = J1\(u1 - ω1 × (J1 * ω1))
    ω2dot = J2\(u2 - ω2 × (J2 * ω2))
    q1dot = Rotations.kinematics(q1, ω1)
    q2dot = Rotations.kinematics(q2, ω2)
    SA[
        q1dot[1], q1dot[2], q1dot[3], q1dot[4],
        ω1dot[1], ω1dot[2], ω1dot[3],
        q2dot[1], q2dot[2], q2dot[3], q2dot[4],
        ω2dot[1], ω2dot[2], ω2dot[3],
    ]
end

RobotDynamics.control_dim(::SatellitePair) = 6

RobotDynamics.LieState(::SatellitePair{R}) where R = RobotDynamics.LieState(R, (0,3,3))

J1 = SMatrix{3,3}(Diagonal(fill(1.0, 3)))
J2 = SMatrix{3,3}(Diagonal(fill(2.0, 3)))
model = SatellitePair{UnitQuaternion{Float64}, Float64}(J1, J2)
x,u = rand(model)
@test norm(x[1:4]) ≈ 1
@test norm(x[8:11]) ≈ 1
s = RobotDynamics.LieState(model)
@test length(s) == 14

ω1 = x[5:7]
ω2 = x[12:14]
q1 = UnitQuaternion(x[1:4],false)
q2 = UnitQuaternion(x[8:11],false)
@test all(RobotDynamics.vec_states(s, x)[2:3] .≈ (ω1, ω2))
@test all(RobotDynamics.rot_states(s, x) .≈ (q1,q2))

model = SatellitePair{MRP{Float64}, Float64}(J1, J2)
x,u = rand(model)
@test length(x) == 12
s = RobotDynamics.LieState(model)
@test length(s) == 12

ω1 = x[4:6]
ω2 = x[10:12]
q1 = MRP(x[1:3]...)
q2 = MRP(x[7:9]...)
@test all(RobotDynamics.vec_states(s, x)[2:3] .≈ (ω1, ω2))
@test all(RobotDynamics.rot_states(s, x) .≈ (q1,q2))

ext: jl


text:
#
println("Ai) Test of the jj-LS transformation of levels from a given multiplet.")
#
wa = Atomic.Computation("jj-LS level transformation",  Nuclear.Model(26.);
                        configs=[Configuration("[Ne] 3s^2 3p^4 3d")],
                        asfSettings=AsfSettings(true, false, "meanDFS", "hydrogenic", "xxx", [1],     0, 1.0e-6, JAC.Subshell[],
                                                true, false, NoneQed(), "yyy", ManyElectron.LSjjSettings(true),
                                                false, [1,2,3,4], false, JAC.LevelSymmetry[] ) )

wb = perform(wa)

##                        configs=[Configuration("[Ne] 3s^2 3p^4"), Configuration("[Ne] 3s^2 3d^4")],  , Configuration("1s^2 2s 2p^6")


text:
using Pkg
Pkg.add("IJulia")

max_stars_repo_name: snigdhayan/Utilities

max_issues_repo_name: snigdhayan/Utilities

max_forks_repo_name: snigdhayan/Utilities


text:
using Plots, Printf, LinearAlgebra

function compute_flux!(qH, H, D, dx, nx)
    Threads.@threads for ix=1:nx
        if (ix<=nx-1)  qH[ix] = -D*(H[ix+1]-H[ix])/dx  end
    end
    return
end

function compute_rate!(ResH, dHdt, H, Hold, qH, dt, damp, dx, nx)
    Threads.@threads for ix=1:nx
        if (2<=ix<=nx-1)  ResH[ix-1] = -(H[ix] - Hold[ix])/dt -(qH[ix]-qH[ix-1])/dx  end
        if (2<=ix<=nx-1)  dHdt[ix-1] = ResH[ix-1] + damp*dHdt[ix-1]  end
    end
    return
end

function compute_update!(H, dHdt, dtau, nx)
    Threads.@threads for ix=1:nx
        if (2<=ix<=nx-1)  H[ix] = H[ix] + dtau*dHdt[ix-1]  end
    end
    return
end

@views function diffusion_1D()
    @show Threads.nthreads()
    # Physics
    lx     = 10.0       # domain size
    D      = 1.0        # diffusion coefficient
    ttot   = 0.6        # total simulation time
    dt     = 0.1        # physical time step
    # Numerics
    nx     = 256        # numerical grid resolution
    tol    = 1e-6       # tolerance
    itMax  = 1e5        # max number of iterations
    damp   = 1-41/nx    # damping (this is a tuning parameter, dependent on e.g. grid resolution)
    # Derived numerics
    dx     = lx/nx      # grid size
    dtau   = (1.0/(dx^2/D/2.1) + 1.0/dt)^-1 # iterative timestep
    xc     = LinRange(dx/2, lx-dx/2, nx)
    # Array allocation
    qH     = zeros(nx-1)
    dHdt   = zeros(nx-2)
    ResH   = zeros(nx-2)
    # Initial condition
    H0     = exp.(-(xc.-lx/2).^2)
    Hold   = copy(H0)
    H      = copy(H0)
    t = 0.0; it = 0; ittot = 0
    # Physical time loop
    while t<ttot
        iter = 0; err = 2*tol
        # Pseudo-transient iteration
        while err>tol && iter<itMax
            compute_flux!(qH, H, D, dx, nx)
            compute_rate!(ResH, dHdt, H, Hold, qH, dt, damp, dx, nx)
            compute_update!(H, dHdt, dtau, nx)
            iter += 1; err = norm(ResH)/length(ResH)
        end
        ittot += iter; it += 1; t += dt
        Hold .= H
    end
    # Analytic solution
    Hana = 1/sqrt(4*(ttot+1/4)) * exp.(-(xc.-lx/2).^2 /(4*(ttot+1/4)))
    @printf("Total time = %1.2f, time steps = %d, iterations tot = %d, error vs analytic = %1.2e \n", round(ttot, sigdigits=2), it, ittot, norm(H-Hana))
    # Visualise
    plot(xc, H0, linewidth=3); display(plot!(xc, H, legend=false, framestyle=:box, linewidth=3, xlabel="lx", ylabel="H", title="damped diffusion (nt=$it, iters=$ittot)"))
    return
end

diffusion_1D()


text:
#
# Custom NDTensors.similar implementation
# More extensive than Base.similar
#

# TODO: define
# Base.similar(t::Tensor, args...) = similar(t, args...)

# In general define NDTensors.similar = Base.similar
similar(a::Array, args...) = Base.similar(a, args...)

# XXX: this is type piracy but why doesn't base have something like this?
# This type piracy is pretty bad, consider making an internal `NDTensors.similar` function.
# Currently Base gives:
# julia> similar(Matrix{Float64}, 2)
# ERROR: MethodError: no method matching Matrix{Float64}(::UndefInitializer, ::Tuple{Int64})
# [...]
similar(::Type{<:Array{T}}, dims) where {T} = Array{T,length(dims)}(undef, dims)

# XXX: this is type piracy but why doesn't base have something like this?
function similar(::Type{ArrayT}, ::Type{ElT}, size) where {ArrayT<:AbstractArray,ElT}
  return similar(similartype(ArrayT, ElT), size)
end

#
# similartype returns the type of the object that would be returned by `similar`
#

# TODO: extend to AbstractVector, returning the same type as `Base.similar` would
# (make sure it handles views, CuArrays, etc. correctly)
# This is to help with code that is generic to different storage types.
similartype(::Type{<:Array{<:Any,N}}, eltype::Type) where {N} = Array{eltype,N}

#similartype(::Type{LinearAlgebra.Adjoint{Float64, Matrix{Float64}}}, ::Type{Float64})

# A union type of AbstractArrays that are wrappers that define `Base.parent`.
#const WrappedArray = Union{ReshapedArray,Adjoint,Transpose,SubArray}
const WrappedArray{T,AW} = Union{
  ReshapedArray{T,<:Any,AW},
  Transpose{T,AW},
  Adjoint{T,AW},
  Symmetric{T,AW},
  Hermitian{T,AW},
  UpperTriangular{T,AW},
  LowerTriangular{T,AW},
  UnitUpperTriangular{T,AW},
  UnitLowerTriangular{T,AW},
  Diagonal{T,AW},
}

parenttype(::Type{<:WrappedArray{<:Any,P}}) where {P} = P

function similartype(::Type{ArrayT}, eltype::Type, dims::Tuple) where {ArrayT<:WrappedArray}
  return similartype(parenttype(ArrayT), eltype, dims)
end

function similartype(::Type{ArrayT}, eltype::Type) where {ArrayT<:WrappedArray}
  return similartype(parenttype(ArrayT), eltype)
end

function similartype(::Type{ArrayT}) where {ArrayT<:WrappedArray}
  return similartype(parenttype(ArrayT))
end


text:
using MatrixMarket
using SuiteSparse
using BenchmarkTools
using MAT
function solveMatrix(x...)

    #Creo cellArray
    matrices = Array{Any}(8)

    #Array per i risultati
    allTimes = []
    allErrors = []
    allSizes = []
    allMems = []

    #Carico le matrici FEMLAB nel cellArray
    print("Loading data...")
    openedMat = matread("../matrices/mat/one20414.mat")
    matrices[1] = openedMat["A"]
    openedMat = matread("../matrices/mat/two367.mat")
    matrices[2] = openedMat["A"]
    openedMat = matread("../matrices/mat/three13514.mat")
    matrices[3] = openedMat["A"]
    openedMat = matread("../matrices/mat/four85623.mat")
    matrices[4] = openedMat["A"]
    openedMat = matread("../matrices/mat/five415.mat")
    matrices[5] = openedMat["A"]
    openedMat = matread("../matrices/mat/six12504.mat")
    matrices[6] = openedMat["A"]
    openedMat = matread("../matrices/mat/seven29067.mat")
    matrices[7] = openedMat["A"]
    openedMat = matread("../matrices/mat/eigth42930.mat")
    matrices[8] = openedMat["A"]

    #Calcola quanti parametri in input ci sono
    howManyOptional = length(x)

    #Matrici in input - se ce ne sono
    if howManyOptional>0
        inputMatrices = Array{Any}(howManyOptional)
        for ll=1 : howManyOptional
            openedMat = matread(x[ll])
            inputMatrices[ll] = openedMat["A"]
        end
    end

    #Codice da eseguire se ci sono matrici in input
    if (howManyOptional>0)
        print("\nProcessing input matrices...")
        (allMems, allTimes, allErrors, allSizes) = interCalc(inputMatrices, allMems, allTimes, allErrors, allSizes)
    else
        print("\nProcessing FEMLAB matrices...")
        (allMems, allTimes, allErrors, allSizes) = interCalc(matrices, allMems, allTimes, allErrors, allSizes)
    end

    #Lunghezza
    if howManyOptional>0
        numDiMatrici = howManyOptional;
    else
        numDiMatrici = 8;
    end

    #Creo matrici contenti i valori dei grafici
    #Ordino anche in modo crescento rispetto alla size
    print("Collecting results...\n")
    firstTime = [allSizes ; allTimes];
    firstTime = reshape(firstTime,numDiMatrici,2);
    firstTimeSor = sortrows(firstTime,by=y->(y[1]))
    firstTimeSor = reshape(firstTimeSor,numDiMatrici,2);

    secondErrors = [allSizes ; allErrors];
    secondErrors = reshape(secondErrors,numDiMatrici,2);
    secondErrorsSor = sortrows(secondErrors,by=y->(y[1]))
    secondErrorsSor = reshape(secondErrorsSor,numDiMatrici,2);

    thirdRAM = [allSizes ; allMems];
    thirdRAM = reshape(thirdRAM,numDiMatrici,2);
    thirdRAMSor = sortrows(thirdRAM,by=y->(y[1]));
    thirdRAMSor = reshape(thirdRAMSor,numDiMatrici,2);

	print(minimum(firstTimeSor[:,2]));
	print("\n");
	print(maximum(firstTimeSor[:,2]));
	print("\n");
	print(sum(firstTimeSor[:,2]));
	print("\n_________________________________");
	print("\n\n\n");

	print(minimum(secondErrorsSor[:,2]));
	print("\n");
	print(maximum(secondErrorsSor[:,2]));
	print("\n");
	print(sum(secondErrorsSor[:,2]));
	print("\n_________________________________");
	print("\n\n\n");
	print(minimum(thirdRAMSor[:,2]));
	print("\n");
	print(maximum(thirdRAMSor[:,2]));
	print("\n");
	print(sum(thirdRAMSor[:,2]));
print("\n_________________________________");
    return firstTimeSor, secondErrorsSor, thirdRAMSor;

end

function interCalc(matrices, allMems, allTimes, allErrors, allSizes)
     for o=1 : size(matrices,1)
        #gc(); workspace();
        tic();

        #Salvo la soluzione esatta xe del sistema
        xe = ones(size(matrices[o],1));

        #Calcolo del right hand side
        b = matrices[o]*xe;
		matrices[o] = lufact(matrices[o]);

        #Risoluzione
        profiling = @timed(matrices[o]\b);
        profiling = profiling[3];
        x = matrices[o]\b;
        print("Processing matrix ",size(matrices[o],1));

        #Errore relativo
        errRel = norm(x-xe)/norm(xe);

        #Salvo dati utili
        elapsed = toc();
        allTimes = [allTimes ; elapsed];
        allSizes = [allSizes ; size(matrices[o],1)];
        allErrors = [allErrors ; errRel];
        #allMems = [allMems];
        allMems = [allMems ; profiling];

    end

    return allMems, allTimes, allErrors, allSizes;

end


text:
using XLNet
using Documenter

DocMeta.setdocmeta!(XLNet, :DocTestSetup, :(using XLNet); recursive=true)

makedocs(;
    modules=[XLNet],
    authors="Sambhaw Kumar, <vampiresambhaw@gmail.com>",
    repo="https://github.com/SambhawDrag/XLNet.jl/blob/{commit}{path}#{line}",
    sitename="XLNet.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://SambhawDrag.github.io/XLNet.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/SambhawDrag/XLNet.jl",
)


text:
# coefficients a_0, ..., a_9
let
  pdTableI = [ -0.5000000002 ;
                0.3333333343 ;
               -0.2499998565 ;
                0.1999997049 ;
               -0.1666848753 ;
                0.1428833286 ;
               -0.1241963125 ;
                0.1101687109 ;
               -0.1142650302 ;
                0.1055093006 ]

  # factorials of 0--9
  pdFactorials = [ 1; 1; 2; 6; 24; 120; 720; 5040; 40320; 362880 ]

  invSqrt2Pi = 1/sqrt(2*π)
  inv12 = 1/12
  inv360 = 1/360
  inv1260 = 1/1260

  @inline function procedureF(K::Int64, μ::Float64, s::Float64, ω::Float64,
    c0::Float64, c1::Float64, c2::Float64, c3::Float64)
    ## 1. Poisson probabilities
    px::Float64 = 0.0
    py::Float64 = 0.0
    if K < 10
      px = -μ
      py = μ^K/pdFactorials[K+1]
    else
      KSq::Float64 = K*K
      KCubed::Float64 = K*KSq
      KFived::Float64 = KCubed * KSq
      δ::Float64 = inv12/K - inv360/KCubed + inv1260/KFived
      V::Float64 = (μ-K)/K
      if abs(V) > 0.25
        px = K*log(1+V) - μ + K - δ
      else
        ws::Float64 = 0.0
        Vj::Float64 = 1.0
        for j = 0:9
          ws += pdTableI[j+1]*Vj
          Vj *= V
        end
        px = K*V*V*ws - δ
      end
      py = 1/sqrt(2*π*K)
    end
    ## 2. Discrete normal probabilities
    X::Float64 = (K - μ + 0.5)/s
    XSq::Float64 = X*X
    fx::Float64 = -0.5*XSq ## typo in paper omits minus sign!
    fy::Float64 = ω*(((c3*XSq + c2)*XSq + c1)*XSq + c0)
    return px, py, fx, fy
  end

  inv24 = 1/24
  inv7 = 1/7

  @inline function caseA(μ::Float64, rng::RNG) where RNG <: AbstractRNG
    s::Float64 = sqrt(μ)
    d::Float64 = 6*μ*μ
    L::Int64 = floor(Int64, μ - 1.1484)
    ## Normal sample
    G::Float64 = μ + s*randn(rng)
    U::Float64 = 0.0
    K::Int64 = 0
    if G >= 0
      K = floor(Int64, G)
      ## Immediate acceptance
      if K >= L return K end
      U = rand(rng)
      ## Squeeze acceptance
      if d*U >= (μ - K)^3 return K end
    end
    ## Preparations for Q and H
    ω::Float64 = invSqrt2Pi/s
    b1::Float64 = inv24/μ
    b2::Float64 = 0.3*b1*b1
    c3::Float64 = inv7*b1*b2
    c2::Float64 = b2 - 15*c3
    c1::Float64 = b1 - 6*b2 + 45*c3
    c0::Float64 = 1 - b1 + 3*b2 - 15*c3
    c::Float64 = 0.1069/μ
    px::Float64 = 0.0
    py::Float64 = 0.0
    fx::Float64 = 0.0
    fy::Float64 = 0.0
    if G >= 0
      px, py, fx, fy = procedureF(K, μ, s, ω, c0, c1, c2, c3)
      ## Quotient acceptance
      if fy*(1-U) <= py*exp(px-fx) return K end
    end
    ## Double exponential sample
    while (true)
      E::Float64 = randexp(rng)
      U = 2*rand(rng) - 1
      T::Float64 = 1.8 + E*sign(U)
      if T <= -0.6744
        continue
      end
      K = floor(Int64, μ + s*T)
      px, py, fx, fy = procedureF(K, μ, s, ω, c0, c1, c2, c3)
      ## Hat acceptance
      if c*abs(U) <= py*exp(px+E) - fy*exp(fx+E)
        return K
      end
    end
  end

  ## modification: no table, just inversion.
  @inline function caseB(μ::Float64, rng::RNG) where RNG <: AbstractRNG
    p::Float64 = exp(-μ)
    q::Float64 = p
    U::Float64 = rand(rng)
    K::Int64 = 0
    while U > q
      K += 1
      p *= μ/K
      q += p
    end
    return K
  end

  ## The algorithm from
  ## Ahrens, J.H. and Dieter, U., 1982. Computer generation of Poisson deviates
  ## from modified normal distributions. ACM Transactions on Mathematical
  ## Software (TOMS), 8(2), pp. 163-179.
  ## This is suitable when μ changes every invocation. Modifications or
  ## alternative algorithms should be used for the case where many Poisson
  ## variates are required for the same μ
  global @inline function samplePoisson(μ::Float64,
      rng::RNG = GLOBAL_RNG) where RNG <: AbstractRNG
    if μ >= 10.0
      return caseA(μ, rng)
    else
      return caseB(μ, rng)
    end
  end
end


text:
#####
##### Codegen/inference integration
#####

function code_cache(mxi::MixtapeInterpreter)
    return WorldView(get_cache(typeof(mxi.inner)), get_world_counter(mxi))
end

function cpu_invalidate(replaced, max_world)
    cache = get_cache(NativeInterpreter)
    invalidate(cache, replaced, max_world, 0)
    return nothing
end

function cpu_cache_lookup(mi, min_world, max_world)
    wvc = WorldView(get_cache(NativeInterpreter), min_world, max_world)
    return Core.Compiler.get(wvc, mi, nothing)
end

function cpu_infer(ctx, mi, min_world, max_world)
    wvc = WorldView(get_cache(NativeInterpreter), min_world, max_world)
    interp = MixtapeInterpreter(ctx, NativeInterpreter(min_world))
    ret = infer(wvc, mi, interp)
    if debug(ctx) && !isempty(interp.errors)
        println("\nEncountered the following non-critical errors during interception:")
        for k in interp.errors
            println(k)
        end
        println("\nThis may imply that your transform failed to operate on some calls.")
        println("\t\t\t  ______________\t\t\n")
    end
    return ret
end

function infer(wvc, mi, interp)
    src = Core.Compiler.typeinf_ext_toplevel(interp, mi)
    ret = Any
    meth = mi.def
    try
        fn = resolve(GlobalRef(mi.def.module, mi.def.name))
        as = map(resolve, mi.specTypes.parameters[2:end])
        ret = Core.Compiler.return_type(interp, fn, as)
    catch e
        push!(interp, e)
    end
    @assert Core.Compiler.haskey(wvc, mi)
    ci = Core.Compiler.getindex(wvc, mi)
    if ci !== nothing && ci.inferred === nothing
        ci.inferred = src
    end
    return ret
end

struct InvokeException <: Exception
    name::Any
    mod::Any
    file::Any
    line::Any
end
function Base.show(io::IO, ie::InvokeException)
    print("@ ($(ie.file), L$(ie.line))\n")
    return print("| (Found call to invoke): $(ie.mod).$(ie.name)\n")
end

function detect_invoke(b, linfo)
    meth = linfo.def
    for (v, st) in b
        st isa Expr || continue
        st.head == :call || continue
        st.args[1] == invoke || continue
        return InvokeException(meth.name, meth.module, meth.file, meth.line)
    end
    return nothing
end

# Replace usage sited of `retrieve_code_info`, OptimizationState is one such, but in all interesting use-cases
# it is derived from an InferenceState. There is a third one in `typeinf_ext` in case the module forbids inference.
function InferenceState(result::InferenceResult, cached::Bool, interp::MixtapeInterpreter)
    src = retrieve_code_info(result.linfo)
    mi = result.linfo
    meth = mi.def
    try
        fn = resolve(GlobalRef(meth.module, meth.name))
        as = map(resolve, result.argtypes[2:end])
        if debug(interp.ctx)
            print("@ ($(meth.file), L$(meth.line))\n")
            print("| beg (inf): $(meth.module).$(fn)\n")
        end
        if allow(interp.ctx, meth.module, fn, as...)
            b = CodeInfoTools.Builder(src, length(result.argtypes[2:end]))
            b = transform(interp.ctx, b)
            e = detect_invoke(b, result.linfo)
            if e != nothing
                push!(interp, e)
            end
            src = finish(b)
        end
    catch e
        push!(interp, e)
    end
    src === nothing && return nothing
    validate_code_in_debug_mode(result.linfo, src, "lowered")
    return InferenceState(result, src, cached, interp)
end

function cpu_compile(ctx, mi, world)
    params = Base.CodegenParams(; track_allocations=false, code_coverage=false,
                                prefer_specsig=true,
                                lookup=@cfunction(cpu_cache_lookup, Any, (Any, UInt, UInt)))

    # generate IR
    # TODO: Instead of extern policy integrate with Orc JIT

    # populate the cache
    if cpu_cache_lookup(mi, world, world) === nothing
        cpu_infer(ctx, mi, world, world)
    end

    # TODO: actual return type.
    rt = Any

    native_code = ccall(:jl_create_native, Ptr{Cvoid},
                        (Vector{Core.MethodInstance}, Base.CodegenParams, Cint), [mi],
                        params, 1) #=extern policy=#
    @assert native_code != C_NULL
    llvm_mod_ref = ccall(:jl_get_llvm_module, LLVM.API.LLVMModuleRef, (Ptr{Cvoid},),
                         native_code)
    @assert llvm_mod_ref != C_NULL
    llvm_mod = LLVM.Module(llvm_mod_ref)

    # get the top-level code
    code = cpu_cache_lookup(mi, world, world)

    # get the top-level function index
    llvm_func_idx = Ref{Int32}(-1)
    llvm_specfunc_idx = Ref{Int32}(-1)
    ccall(:jl_get_function_id, Nothing, (Ptr{Cvoid}, Any, Ptr{Int32}, Ptr{Int32}),
          native_code, code, llvm_func_idx, llvm_specfunc_idx)
    @assert llvm_func_idx[] != -1
    @assert llvm_specfunc_idx[] != -1

    # get the top-level function
    llvm_func_ref = ccall(:jl_get_llvm_function, LLVM.API.LLVMValueRef,
                          (Ptr{Cvoid}, UInt32), native_code, llvm_func_idx[] - 1)
    @assert llvm_func_ref != C_NULL
    llvm_func = LLVM.Function(llvm_func_ref)
    llvm_specfunc_ref = ccall(:jl_get_llvm_function, LLVM.API.LLVMValueRef,
                              (Ptr{Cvoid}, UInt32), native_code, llvm_specfunc_idx[] - 1)
    @assert llvm_specfunc_ref != C_NULL
    llvm_specfunc = LLVM.Function(llvm_specfunc_ref)

    return (rt, llvm_specfunc, llvm_func, llvm_mod)
end

function method_instance(@nospecialize(f), @nospecialize(tt), world)
    # get the method instance
    meth = which(f, tt)
    sig = Base.signature_type(f, tt)::Type
    (ti, env) = ccall(:jl_type_intersection_with_env, Any, (Any, Any), sig,
                      meth.sig)::Core.SimpleVector
    meth = Base.func_for_method_checked(meth, ti, env)
    return ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance},
                 (Any, Any, Any, UInt), meth, ti, env, world)
end

function codegen(ctx::CompilationContext, @nospecialize(f), @nospecialize(tt);
                 world=Base.get_world_counter())
    return cpu_compile(ctx, method_instance(f, tt, world), world)
end


text:
struct VanDerWaals{S,T,A} <: CubicModel
    type::S
    tc::T
    pc::T
    ω::T
    mw::T
    vc::Union{T,Nothing}
    _a::T
    _b::T
    aij::A

    function VanDerWaals(tc,pc,ω,mw,vc=nothing,aij = nothing)
        if length(tc) == 1
            type = SINGLE()

        else
            type = MULTI()
        end
        S = typeof(type)
        tc,pc,ω,mw = promote(tc,pc,ω,mw)
        T = typeof(tc)
        A = typeof(aij)
        _a = 0.421875 .*((RGAS .*tc).^2) ./pc
        _b = 0.125 .*((RGAS .*tc)) ./pc
        return new{S,T,A}(type,tc,pc,ω,mw,vc,_a,_b,aij)
    end
end

volume_solver_type(model::VanDerWaals) = CubicRoots()
single_sat_Approx(model::VanDerWaals{SINGLE}) =VdWSatApprox(model)
mol_density(model::VanDerWaals{SINGLE},::CriticalPoint,unit=u"mol/(m^3)") = convert_unit(u"mol/L",unit,inv(only(model.vc)))
pressure(model::VanDerWaals{SINGLE},::CriticalPoint,unit=u"Pa") = convert_unit(u"Pa",unit,only(model.pc))
temperature(model::VanDerWaals{SINGLE},::CriticalPoint,unit=u"K") = convert_unit(u"K",unit,only(model.tc))
mol_volume(model::VanDerWaals{SINGLE},::CriticalPoint,unit=u"m^3/mol") = convert_unit(u"m^3/mol",unit,only(model.vc))
acentric_factor(model::VanDerWaals{SINGLE}) = only(model.ω)
molecular_weight(model::VanDerWaals{SINGLE}) = only(model.mw)
molecular_weight(model::VanDerWaals{MULTI}) = model.mw

function mol_density(model::VanDerWaals{MULTI},::CriticalPoint,unit=u"mol/(m^3)")
    return convert_unit.(u"mol/L",unit,1 ./ model.vc)
end

function pressure(model::VanDerWaals{MULTI},::CriticalPoint,unit=u"Pa")
    return convert_unit.(u"Pa",unit,model.pc)
end

function temperature(model::VanDerWaals{MULTI},::CriticalPoint,unit=u"K")
    return convert_unit.(u"K",unit,model.tc)
end

function mol_volume(model::VanDerWaals{MULTI},::CriticalPoint,unit=u"m^3/mol")
    return convert_unit.(u"m^3/mol",unit,model.vc)
end

function acentric_factor(model::VanDerWaals{MULTI})
    return model.ω
end

function VanDerWaals(;tc,pc,ω,mw,vc=nothing,aij=nothing)
    return VanDerWaals(tc,pc,ω,mw,vc,aij)
end

function VanDerWaals(model::ThermoModel)
    tc = temperature(model,CriticalPoint())
    pc = pressure(model,CriticalPoint())
    ω = acentric_factor(model)
    mw = molecular_weight(model)
    return VanDerWaals(tc,pc,ω,mw)
end

function cubic_ab(mt::SinglePT,model::VanDerWaals{SINGLE},v,t)
    a = only(model._a)
    b = only(model._b)
    return a,b
end

function cubic_ab(mt::MultiPT,model::VanDerWaals{MULTI},p,t,x)

    #two options to introduce alpha:
    #here: it will allocate, but less ops
    ai = model._a
    bi = model._b
    b = dot(bi,x)
    #here: it will not allocate, but more ops
    mixrule= (ai,aj)->sqrt(ai*aj)
    a = mixing_rule(mixrule, x, ai,model.aij)
    return a,b
end

function cubic_abp(mt::SingleVT,model::VanDerWaals{SINGLE},v,t)
    a,b = cubic_ab(QuickStates.pt(),model,v,t) #v is ignored
    p = RGAS*t/(v-b) - a/(v*v)
    return a,b,p
end

function cubic_abp(mt::MultiVT,model::VanDerWaals{MULTI},v,t,x)
    a,b = cubic_ab(QuickStates.ptx(),model,v,t,x) #v is ignored
    p =  RGAS*t/(v-b) - a/(v*v)
    return a,b,p
end

function fugacity_coeff_impl(mt::SingleVT,model::VanDerWaals{SINGLE},v,t)
    a,b,p =  cubic_abp(mt,model,v,t)
    RTinv = 1/(RGAS*t)
    A = a*p*RTinv*RTinv
    B = b*p*RTinv
    z = p*v*RTinv
     _1 = one(z)
     logϕ = z - _1 - log(z-B) - A/z
end

function  αR_impl(mt::MultiVT,model::VanDerWaals{MULTI},rho,t,x)
    R = RGAS
    RTinv = 1/(RGAS*t)
    v = inv(rho)
    a,b,p =  cubic_abp(mt,model,v,t,x)
    -log(1-b*rho) - a*rho*RTinv
end

function cubic_poly(mt::SinglePT,model::VanDerWaals{SINGLE},p,t)
    a,b = cubic_ab(QuickStates.pt(),model,p,t)
    RTinv = 1/(RGAS*t)
    A = a*p*RTinv*RTinv
    B = b*p*RTinv
    _1 = one(a)
    return (-A*B, A, -B-_1, _1)
end

function cubic_poly(mt::MultiPT,model::VanDerWaals{MULTI},p,t,x)
    a,b = cubic_ab(QuickStates.ptx(),model,p,t,x)
    RTinv = 1/(RGAS*t)
    A = a*p*RTinv*RTinv
    B = b*p*RTinv
    _1 = one(a)
    return (-A*B, A, -B-_1, _1)
end

struct VdWSatApprox{M} <: SaturationModel
    model::M
    function VdWSatApprox(model::CubicModel)
        T = typeof(model)
        return new{T}(model)
    end
end

function VdWSatApprox(tc,pc,mw=copy(pc),vc=nothing,ω=nothing,aij = nothing)
    model = VanDerWaals(tc,pc,mw,vc,ω,aij)
    return VdWSatApprox(model)
end

function VdWSatApprox(;tc,pc,mw=copy(pc),vc=nothing,ω=nothing,aij=nothing)
    model =  VanDerWaals(tc,pc,mw,vc,ω,aij)
    return VdWSatApprox(model)
end

function pressure_impl(mt::SingleSatT,model::VdWSatApprox,t)
    a = 4.406664258927600
    b = 2.205610041969020
    c = 0.243757663628277
    d = -0.206185849671953
    e = -4.650419726136550
    f = 0.019582516700758
    tc = only(model.model.tc)
    pc = only(model.model.pc)
    tr = t/tc
    atr = tc/t
    num = evalpoly(atr,(a,c,e))
    denom = evalpoly(atr,(1.0,b,d,f))
    pr = exp(num/denom)*tr
    p = pr*pc
end

initial_temperature(model::VdWSatApprox,p) = critical_sat_interpolation(model,p)


text:
### Visualization routines to aid development
# This entire implementation is extremely hacky...
# Ideally this should be eventually replaced by something based on RecipesBase
#
# Written by Boris Bonev, Nov. 2020

## routine for plotting ranks
function plotranks(hssA::HssMatrix; cutoff_level=3)
  m,n = size(hssA)
  aspect = m/n
  plot(yflip=true, showaxis=false, size = (400, 400*aspect))
  xticks = [1]; yticks = [1]
  yticks, xticks = _plotranks!(hssA,1,1,xticks,yticks,1,cutoff_level)
  append!(yticks, n+1)
  append!(xticks, m+1)
  plot!(aspect_ratio=:equal)
  plot_ref = plot!(xticks=xticks, yticks=yticks, xmirror=true)
  return plot_ref
end

# auxiliary routine for drawing a rectangle
_rectangle(w, h, x, y) = Shape([0,w,w,0] .+ x , [0,0,h,h] .+ y)

# recursive plotting routine
function _plotranks!(hssA::HssMatrix, co, ro, cticks, rticks, level, cl)
  if isleaf(hssA)
    m, n = size(hssA)
    plot!(_rectangle(n,m,co,ro), color=:tomato, label=false)
  else
    m1, n1 = hssA.sz1
    m2, n2 = hssA.sz2
    # plot diagonal blocks and ticks
    _plotranks!(hssA.A11, co, ro, cticks, rticks, level+1, cl)
    if level < cl
      append!(rticks, ro+m1)
      append!(cticks, co+n1)
    end
    _plotranks!(hssA.A22, co+n1, ro+m1, cticks, rticks, level+1, cl)
    # plot off-diagonal blocks
    plot!(_rectangle(n2, m1, co+n1, ro), color=:aliceblue, label=false)
    if level ≤ cl annotate!((co+n1+0.5*n2, ro+0.5*m1, text(max(size(hssA.B12)...), 8))) end
    plot!(_rectangle(n1, m2, co, ro+m1), color=:aliceblue, label=false)
    if level ≤ cl annotate!((co+0.5*n1, ro+m1+0.5*m2, text(max(size(hssA.B21)...), 8))) end
  end
  return rticks, cticks
end

# convenience function imitating the behavior of pcolor in matlab
pcolor(A; args...) = heatmap(A, xlims=(1,size(A,2)), ylims=(1,size(A,1)), yflip=true, aspect_ratio=:equal; args...);


text:
using RecipesBase
#using Plots
#gr()

using Mapper2.MapperGraphs
#pyplot()

################################################################################
# Post routing plotting.
################################################################################
struct Polygon
    x::Vector{Float64}
    y::Vector{Float64}
    fill::Symbol
end

# ------- #
# Methods #
# ------- #

# Given a box, give X and Y vectors to trace the outline of the box.
getx(d::Polygon) = d.x
gety(d::Polygon) = d.y

# Encoding for a route through the architecture.
struct DrawRoute
    x       ::Vector{Float64}
    y       ::Vector{Float64}
    color   ::Symbol
end

@recipe function f(map::Map; plot_route = true)
    # Unpack map
    #map = r.map
    #plot_route = r.plot_route

    # Set up parameters
    spacing = 10
    tilesize = 20

    # Build boxes and routes
    boxes = getboxes(map)
    routes = plot_route ? getroutes(map) : getlines(map)

    # Set up plot attributes
    legend := false
    ticks  := nothing
    grid   := false
    yflip  := true

    # Plot boxes
    seriestype := :shape

    for box in boxes
        @series begin
            #linecolor := :black
            c := box.fill

            # Get x,y coordinates from box
            x = getx(box)
            y = gety(box)
            x, y
        end
    end

    seriestype := :path
    linewidth  := 2

    # Plot routes
    for route in routes
        @series begin
            linecolor := route.color
            x = route.x
            y = route.y
            x,y
        end
    end
end

function getboxes(m::Map{2})
    a = m.toplevel
    style = a.metadata["style"]

    # Create draw boxes for each tile in the array.
    boxes = Polygon[]
    for (name, child) in a.children
        addr = getaddress(a, name)
        # scale x,y
        if haskey(child.metadata, "shadow_offset")
            addrs = [addr + o for o in child.metadata["shadow_offset"]]
            push!(addrs, addr)
        else
            addrs = [addr]
        end

        # Unpack tuple after manipulation
        y, x = cartesian(style, dim_min(addrs))

        # fill with cyan if box address is used.
        if MapperCore.isused(m, addr)
            fill = :cyan
        else
            fill = :white
        end
        push!(boxes, Polygon(polygon(style, x, y)..., fill))
    end
    return boxes
end

function iscongested(edge::Mapper2.TaskgraphEdge)
    metadata = edge.metadata
    return get(metadata, "iscongested", false)
end

function getroutes(m::Map{2})
    a = m.toplevel
    style = a.metadata["style"]

    routes = DrawRoute[]
    for (edge_index, graph) in enumerate(m.mapping.edges)
        x = Float64[]
        y = Float64[]
        for path in Mapper2.MapperGraphs.linearize(graph)
            # Only look at global port paths.
            isglobalport(path) || continue
            # Get the address from the path.
            address = getaddress(a, path)

            # Create offsets for smooth paths
            metadata = a[path].metadata
            y0, x0 = cartesian(style, address)

            push!(x, x0 + get(metadata, "x", 0.5))
            push!(y, y0 + get(metadata, "y", 0.5))
        end
        # Choose color based on length of path
        edge = getedge(m.taskgraph, edge_index)

        if iscongested(edge)
            color = :red
        elseif length(x) <= 2
            color = :black
        elseif length(x) <= 5
            color = :blue
        else
            color = :green
        end
        # Add this route to the routes vector
        push!(routes, DrawRoute(x, y, color))
    end
    return routes
end

function getlines(m::Map{2})
    a = m.toplevel
    lines = DrawRoute[]
    for edge in getedges(m.taskgraph)

        x = Float64[]
        y = Float64[]
        source = first(getsources(edge))
        dest   = first(getsinks(edge))

        for node in (source, dest)
            # Get the address from the path.
            path = Mapper2.MapperCore.getpath(m.mapping, node)
            address = getaddress(a, path)
            # Create offsets for smooth paths

            # Big offset for macro location in the whole array
            x_offset_big = address[2]*(spacing + tilesize)
            y_offset_big = address[1]*(spacing + tilesize)
            # Small offset for offset within a tile
            x_offset_small = get(a[path].metadata, "x", 0.5) * tilesize
            y_offset_small = get(a[path].metadata, "y", 0.5) * tilesize

            push!(x, x_offset_big + x_offset_small)
            push!(y, y_offset_big + y_offset_small)
        end
        # Default to black color
        color = :black
        # Add this route to the routes vector
        push!(lines, DrawRoute(x,y,color))
    end
    return lines
end


text:
# offsetting polygons

Using Luxor

@png begin
    background("aliceblue")
    p = star(O, 45, 5, 0.5, 0, vertices=true)
    setline(5)
    sethue("midnightblue")
    setdash("dotdotdashed")
    Luxor.poly(p, :stroke, close=true)
    setdash("solid")
    sethue("deeppink4")

    setline(2)
    Luxor.poly(offsetpoly(p, 40), :stroke, close=true)
    Luxor.poly(offsetpoly(p, 45), :stroke, close=true)
    Luxor.poly(offsetpoly(p, 50), :stroke, close=true)
    Luxor.poly(offsetpoly(p, 55), :stroke, close=true)

    sethue("darkcyan")

    Luxor.poly(offsetpoly(p, 30), :stroke, close=true)
    Luxor.poly(offsetpoly(p, 35), :stroke, close=true)
    Luxor.poly(offsetpoly(p, 40), :stroke, close=true)
    Luxor.poly(offsetpoly(p, 45), :stroke, close=true)

    sethue("tan")

    setline(1)
    Luxor.poly(offsetpoly(p, 30), :stroke, close=true)
    Luxor.poly(offsetpoly(p, 25), :stroke, close=true)
    Luxor.poly(offsetpoly(p, 30), :stroke, close=true)
    Luxor.poly(offsetpoly(p, 35), :stroke, close=true)

    sethue("deepskyblue4")

    Luxor.poly(offsetpoly(p, -10), :stroke, close=true)
    Luxor.poly(offsetpoly(p, -15), :stroke, close=true)
    Luxor.poly(offsetpoly(p, -20), :stroke, close=true)
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
insts = OrderedDict()
insts["minimal"] = [
    ((4, true),),
    ((4, false), :SOCExpPSD),
    ]
insts["fast"] = [
    ((6, true),),
    ((6, false), :SOCExpPSD),
    ((8, true),),
    ]
insts["various"] = vcat(insts["fast"], [
    ((8, false), :SOCExpPSD),
    ((10, true),),
    ((10, false), :SOCExpPSD),
    ])
return (RegionOfAttrJuMP, insts)


text:
using EchelleCCFs
using RvSpectMLBase
using Test

@testset "EchelleCCFs.jl" begin

    include("util.jl")

    include("linelist.jl")

    include("mask_shapes.jl")

    include("ccf_plan.jl")

    include("calc_rv.jl")

    include("convenience.jl")

end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
@reexport module RandomWalk

export RandomWalkEnv

using ReinforcementLearningBase

mutable struct RandomWalkEnv <: AbstractEnv
    N::Int
    start::Int
    state::Int
    actions::Vector{Int}
    leftreward::Float64
    rightreward::Float64
    observation_space::DiscreteSpace
    action_space::DiscreteSpace

    RandomWalkEnv(
        ;
        N = 7,
        actions = [-1, 1],
        start = div(N + 1, 2),
        leftreward = -1.0,
        rightreward = 1.0,
    ) =
        new(
            N,
            start,
            start,
            actions,
            leftreward,
            rightreward,
            DiscreteSpace(N),
            DiscreteSpace(length(actions)),
        )
end

RLBase.get_observation_space(env::RandomWalkEnv) = env.observation_space
RLBase.get_action_space(env::RandomWalkEnv) = env.action_space

function (env::RandomWalkEnv)(a::Int)
    env.state = min(max(env.state + env.actions[a], 1), env.N)
    nothing
end

function RLBase.reset!(env::RandomWalkEnv)
    env.state = env.start
    nothing
end

RLBase.observe(env::RandomWalkEnv) =
    (
        state = env.state,
        terminal = env.state == env.N || env.state == 1,
        reward = env.state == env.N ? env.rightreward :
                 (env.state == 1 ? env.leftreward : 0.0),
    )

end


text:
# notebook will be generated
x = 1 + 2


text:
"""
Circular Restricted Three-body Model tests.
"""
module PropagationTests

using LinearAlgebra, Unitful, UnitfulAstro, GeneralAstrodynamics, Test
using DifferentialEquations

@testset verbose=false "R2BP Propagation" begin

    orbit = let planet = Mars
        e = 0.4
        a = 10_000
        i = 0
        Ω = 0
        ω = 0
        ν = 0
        state = KeplerianState(e, a, i, Ω, ω, ν)
        state = CartesianState(cartesian(state, massparameter(planet))...)
        Orbit(state, planet)
    end

    traj = propagate(orbit, period(orbit); save_everystep=false)

    @test traj[end] ≈ state(orbit)

end

@testset verbose=false "CR3BP Propagation" begin

    @test isperiodic(halo(SunEarth;   L=2, Az=0.0005)...)
    @test isperiodic(halo(SunMars;    L=1, Az=0.0005)...)
    @test isperiodic(halo(EarthMoon;  L=2, Az=0.0005)...)
    @test isperiodic(halo(SunJupiter; L=1, Az=0.0005)...)

    orbit, T = halo(SunJupiter; L=2, Az=0)
    @test isapprox.(
        monodromy(orbit, T), [
            1177.6158450389235 -43.31366732210663 0.0 247.08544503455505 227.41064838834146 0.0;
           -1085.995406851377 40.86470821207657 0.0 -227.41064838735218 -209.97647807144125 0.0;
            0.0 0.0 0.9966422601737439 0.0 0.0 -0.09182654535515093;
            3446.1566018075323 -126.52938312132042 0.0 722.7945482654276 666.0424507143235 0.0;
           -2146.972890519174 79.03452084349573 0.0 -450.8572227441975 -413.95658856183604 0.0;
            0.0 0.0 0.07300944628976104 0.0 0.0 0.99664226017845
        ]; atol=1e-2
    ) |> all

end

end # module


text:
using PSCOPF

using .PSCOPFFixtures

using Test
using JuMP
using Dates
using DataStructures
using Printf

@testset verbose=true "test_eod_assessment_call" begin

    function create_instance(limit_1, impositions_1, impositions_2,
                            limit::Float64=35.,
                            logs=nothing)
        network = PSCOPFFixtures.network_2buses(limit=limit)
        # Limitables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "wind_1_1", PSCOPF.Networks.LIMITABLE,
                                                0., 200.,
                                                0., 1.,
                                                Dates.Second(0), Dates.Second(0))
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                0., 200.,
                                                0., 10.,
                                                Dates.Second(0), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_2", "prod_2_1", PSCOPF.Networks.PILOTABLE,
                                                0., 200.,
                                                0., 50.,
                                                Dates.Second(0), Dates.Second(0))

        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()

        #Fixme to args
        assessment_uncertainties = SortedDict("bus_1" => (30, 50),
                                            "bus_2" => (20, 70),
                                            "wind_1_1" => (50, 90)
                                            )

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                    SortedDict{String,PSCOPF.GeneratorState}(), #gen_initial_state
                                    uncertainties, assessment_uncertainties, logs)

        PSCOPF.set_imposition_definitive_value!(PSCOPF.get_tso_actions(context), "prod_1_1", TS[1], impositions_1[1], impositions_1[2])
        PSCOPF.set_imposition_definitive_value!(PSCOPF.get_tso_actions(context), "prod_2_1", TS[1], impositions_2[1], impositions_2[2])
        PSCOPF.set_limitation_value!(PSCOPF.get_tso_actions(context), "wind_1_1", TS[1], limit_1)

        return context
    end

    TS = [DateTime("2015-01-01T11:00:00")]
    ech = DateTime("2015-01-01T07:00:00")
    next_ech = DateTime("2015-01-01T07:30:00")

    #=
                    bus 1                   bus 2
                            |                      |
        (limitable) wind_1_1|       "1_2"          |
        Pmin=0, Pmax=200    |                      |
        Csta=0, Cprop=1     |                      |
        [50-90]             |----------------------|
                            |         35           |
                            |                      |
        (pilotable) prod_1_1|                      |(pilotable) prod_2_1
        Pmin=20, Pmax=200   |                      | Pmin=20, Pmax=200
        Csta=10k, Cprop=10  |                      | Csta=50k, Cprop=50
        [0-0]               |                      | [20-200]
                            |                      |
               load(bus_1)  |                      |load(bus_2)
        [30-50]             |                      | [20-70]
    =#
    @testset verbose=true "test_no_violation" begin
        limit_1 = 75.
        impositions_1 = (0., 0.)
        impositions_2 = (20., 200.)

        context = create_instance(limit_1, impositions_1, impositions_2,
                            35.,)

        assessment = PSCOPF.EODAssessment()
        result = PSCOPF.run(assessment, ech, TS, context)

        @test PSCOPF.is_validated(result)
    end

    #=
                    bus 1                   bus 2
                            |                      |
        (limitable) wind_1_1|       "1_2"          |
        Pmin=0, Pmax=200    |                      |
        Csta=0, Cprop=1     |                      |
        [50-90]             |----------------------|
                            |         35           |
                            |                      |
        (pilotable) prod_1_1|                      |(pilotable) prod_2_1
        Pmin=20, Pmax=200   |                      | Pmin=20, Pmax=200
        Csta=10k, Cprop=10  |                      | Csta=50k, Cprop=50
        [0-0]               |                      | [20-60]
                            |                      |
               load(bus_1)  |                      |load(bus_2)
        [30-50]             |                      | [20-70]
    =#
    @testset verbose=true "test_violation" begin
        limit_1 = 75.
        impositions_1 = (0., 0.)
        impositions_2 = (20., 60.)

        context = create_instance(limit_1, impositions_1, impositions_2,
                            35.,)

        assessment = PSCOPF.EODAssessment()
        result = PSCOPF.run(assessment, ech, TS, context)

        @test !PSCOPF.is_validated(result)
    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
function FieldCartesianOrbit(arg0::FieldOrbit)
    return FieldCartesianOrbit((FieldOrbit,), arg0)
end

function FieldCartesianOrbit(arg0::FieldPVCoordinates, arg1::Frame, arg2::FieldAbsoluteDate, arg3::RealFieldElement)
    return FieldCartesianOrbit((FieldPVCoordinates, Frame, FieldAbsoluteDate, RealFieldElement), arg0, arg1, arg2, arg3)
end

function FieldCartesianOrbit(arg0::TimeStampedFieldPVCoordinates, arg1::Frame, arg2::RealFieldElement)
    return FieldCartesianOrbit((TimeStampedFieldPVCoordinates, Frame, RealFieldElement), arg0, arg1, arg2)
end

function add_kepler_contribution(obj::FieldCartesianOrbit, arg0::PositionAngle, arg1::RealFieldElement, arg2::Vector{RealFieldElement})
    return jcall(obj, "addKeplerContribution", void, (PositionAngle, RealFieldElement, Vector{RealFieldElement}), arg0, arg1, arg2)
end

function get_a(obj::FieldCartesianOrbit)
    return jcall(obj, "getA", RealFieldElement, ())
end

function get_a_dot(obj::FieldCartesianOrbit)
    return jcall(obj, "getADot", RealFieldElement, ())
end

function get_date(obj::FieldOrbit)
    return jcall(obj, "getDate", FieldAbsoluteDate, ())
end

function get_e(obj::FieldCartesianOrbit)
    return jcall(obj, "getE", RealFieldElement, ())
end

function get_e_dot(obj::FieldCartesianOrbit)
    return jcall(obj, "getEDot", RealFieldElement, ())
end

function get_equinoctial_ex(obj::FieldCartesianOrbit)
    return jcall(obj, "getEquinoctialEx", RealFieldElement, ())
end

function get_equinoctial_ex_dot(obj::FieldCartesianOrbit)
    return jcall(obj, "getEquinoctialExDot", RealFieldElement, ())
end

function get_equinoctial_ey(obj::FieldCartesianOrbit)
    return jcall(obj, "getEquinoctialEy", RealFieldElement, ())
end

function get_equinoctial_ey_dot(obj::FieldCartesianOrbit)
    return jcall(obj, "getEquinoctialEyDot", RealFieldElement, ())
end

function get_frame(obj::FieldOrbit)
    return jcall(obj, "getFrame", Frame, ())
end

function get_hx(obj::FieldCartesianOrbit)
    return jcall(obj, "getHx", RealFieldElement, ())
end

function get_hx_dot(obj::FieldCartesianOrbit)
    return jcall(obj, "getHxDot", RealFieldElement, ())
end

function get_hy(obj::FieldCartesianOrbit)
    return jcall(obj, "getHy", RealFieldElement, ())
end

function get_hy_dot(obj::FieldCartesianOrbit)
    return jcall(obj, "getHyDot", RealFieldElement, ())
end

function get_i(obj::FieldCartesianOrbit)
    return jcall(obj, "getI", RealFieldElement, ())
end

function get_i_dot(obj::FieldCartesianOrbit)
    return jcall(obj, "getIDot", RealFieldElement, ())
end

function get_jacobian_wrt_cartesian(obj::FieldOrbit, arg0::PositionAngle, arg1::Vector{Vector{RealFieldElement}})
    return jcall(obj, "getJacobianWrtCartesian", void, (PositionAngle, Vector{Vector{RealFieldElement}}), arg0, arg1)
end

function get_jacobian_wrt_parameters(obj::FieldOrbit, arg0::PositionAngle, arg1::Vector{Vector{RealFieldElement}})
    return jcall(obj, "getJacobianWrtParameters", void, (PositionAngle, Vector{Vector{RealFieldElement}}), arg0, arg1)
end

function get_keplerian_mean_motion(obj::FieldOrbit)
    return jcall(obj, "getKeplerianMeanMotion", RealFieldElement, ())
end

function get_keplerian_period(obj::FieldOrbit)
    return jcall(obj, "getKeplerianPeriod", RealFieldElement, ())
end

function get_le(obj::FieldCartesianOrbit)
    return jcall(obj, "getLE", RealFieldElement, ())
end

function get_le_dot(obj::FieldCartesianOrbit)
    return jcall(obj, "getLEDot", RealFieldElement, ())
end

function get_lm(obj::FieldCartesianOrbit)
    return jcall(obj, "getLM", RealFieldElement, ())
end

function get_lm_dot(obj::FieldCartesianOrbit)
    return jcall(obj, "getLMDot", RealFieldElement, ())
end

function get_lv(obj::FieldCartesianOrbit)
    return jcall(obj, "getLv", RealFieldElement, ())
end

function get_lv_dot(obj::FieldCartesianOrbit)
    return jcall(obj, "getLvDot", RealFieldElement, ())
end

function get_mu(obj::FieldOrbit)
    return jcall(obj, "getMu", RealFieldElement, ())
end

function get_pv_coordinates(obj::FieldOrbit)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, ())
end

function get_pv_coordinates(obj::FieldOrbit, arg0::FieldAbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (FieldAbsoluteDate, Frame), arg0, arg1)
end

function get_pv_coordinates(obj::FieldOrbit, arg0::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (Frame,), arg0)
end

function get_type(obj::FieldCartesianOrbit)
    return jcall(obj, "getType", OrbitType, ())
end

function has_derivatives(obj::FieldCartesianOrbit)
    return jcall(obj, "hasDerivatives", jboolean, ())
end

function interpolate(obj::FieldCartesianOrbit, arg0::FieldAbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", FieldCartesianOrbit, (FieldAbsoluteDate, Stream), arg0, arg1)
end

function interpolate(obj::FieldTimeInterpolable, arg0::FieldAbsoluteDate, arg1::Collection)
    return jcall(obj, "interpolate", FieldTimeInterpolable, (FieldAbsoluteDate, Collection), arg0, arg1)
end

function shifted_by(obj::FieldCartesianOrbit, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldCartesianOrbit, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldCartesianOrbit, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldCartesianOrbit, (jdouble,), arg0)
end

function to_orbit(obj::FieldCartesianOrbit)
    return jcall(obj, "toOrbit", CartesianOrbit, ())
end

function to_string(obj::FieldCartesianOrbit)
    return jcall(obj, "toString", JString, ())
end


text:
Random.seed!(12345)

@testset "Real Linear Solver" begin
  A = rand(3,2);
  x = rand(2);
  b = A*x;

  solvers = linearSolverListReal()

  for solver in solvers
    solverInfo = SolverInfo(Float64)
    S = createLinearSolver(solver, A, iterations=100, solverInfo=solverInfo, shape=(2,1))
    x_approx = solve(S,b)
    @info "Testing solver $solver ...: $x  == $x_approx"
    @test norm(x - x_approx) / norm(x) ≈ 0 atol=0.1
  end

  for solver in solvers
    solverInfo = SolverInfo(Float32)
    S = createLinearSolver(solver, A, iterations=100, solverInfo=solverInfo, shape=(2,1))
    x_approx = solve(S,b)
    @info "Testing solver $solver ...: $x  == $x_approx"
    @test norm(x - x_approx) / norm(x) ≈ 0 atol=0.1
  end
end

@testset "Complex Linear Solver" begin
  A = rand(3,2)+im*rand(3,2);
  x = rand(2)+im*rand(2);
  b = A*x;

  solvers = linearSolverList()

  for solver in solvers
    solverInfo = SolverInfo(ComplexF64)
    S = createLinearSolver(solver, A, iterations=100, solverInfo=solverInfo)
    x_approx = solve(S,b)
    @info "Testing solver $solver ...: $x  == $x_approx"
    @test norm(x - x_approx) / norm(x) ≈ 0 atol=0.1
  end

  for solver in solvers
    solverInfo = SolverInfo(ComplexF32)
    S = createLinearSolver(solver, A, iterations=100, solverInfo=solverInfo)
    x_approx = solve(S,b)
    @info "Testing solver $solver ...: $x  == $x_approx"
    @test norm(x - x_approx) / norm(x) ≈ 0 atol=0.1
  end
end

@testset "General Convex Solver" begin
  # fully sampled operator, image and data
  N = 256
  numPeaks = 5
  F = [1 / sqrt(N)*exp(-2im * π *j*k/N) for j=0:N-1, k=0:N-1]
  x = zeros(N)
  for i = 1:3
    x[rand(1:N)] = rand()
  end
  b = 1 / sqrt(N)*fft(x)

  # random undersampling
  idx = sort(unique(rand(1:N, div(N,2))))
  b = b[idx]
  F = F[idx,:]

  for solver in ["fista","admm"]
    reg = Regularization("L1",1e-3)
    solverInfo = SolverInfo(ComplexF64)
    S = createLinearSolver(solver,F; reg=reg, iterations=200, solverInfo=solverInfo, normalizeReg=false)
    x_approx = solve(S, b)
    @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
    @test x ≈ x_approx rtol=0.1

    reg.λ *= length(b)/norm(b,1)
    scale_F = 1e3 # test invariance to the maximum eigenvalue
    S = createLinearSolver(solver, F .* scale_F; reg=reg, iterations=200, solverInfo=solverInfo, normalizeReg=true)
    x_approx = solve(S, b)
    x_approx .*= scale_F
    @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
    @test x ≈ x_approx rtol=0.1
  end

  for solver in ["fista","admm"]
    reg = Regularization("L1",1e-3)
    solverInfo = SolverInfo(ComplexF32)
    S = createLinearSolver(solver,F; reg=reg, iterations=200, solverInfo=solverInfo, normalizeReg=false)
    x_approx = solve(S, b)
    @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
    @test x ≈ x_approx rtol=0.1

    reg.λ *= length(b)/norm(b,1)
    scale_F = 1e3 # test invariance to the maximum eigenvalue
    S = createLinearSolver(solver, F .* scale_F; reg=reg, iterations=200, solverInfo=solverInfo, normalizeReg=true)
    x_approx = solve(S, b)
    x_approx .*= scale_F
    @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
    @test x ≈ x_approx rtol=0.1
  end

  # test ADMM with option vary_ρ
  solver = "admm"
  reg = Regularization("L1",1.e-3)
  solverInfo = SolverInfo(ComplexF64)
  S = createLinearSolver(solver,F; reg=reg, iterations=200, solverInfo=solverInfo, normalizeReg=false, ρ=1e6, vary_ρ=true, verbose=false)
  x_approx = solve(S, b)
  @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
  @test x ≈ x_approx rtol=0.1

  S = createLinearSolver(solver,F; reg=reg, iterations=200, solverInfo=solverInfo, normalizeReg=false, ρ=1e-6, vary_ρ=true, verbose=false)
  x_approx = solve(S, b)
  @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
  @test x ≈ x_approx rtol=0.1

  # test ADMM with option vary_ρ
  solver = "admm"
  reg = Regularization("L1",1.e-3)
  solverInfo = SolverInfo(ComplexF32)
  S = createLinearSolver(solver,F; reg=reg, iterations=200, solverInfo=solverInfo, normalizeReg=false, ρ=1e6, vary_ρ=true, verbose=false)
  x_approx = solve(S, b)
  @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
  @test x ≈ x_approx rtol=0.1

  S = createLinearSolver(solver,F; reg=reg, iterations=200, solverInfo=solverInfo, normalizeReg=false, ρ=1e-6, vary_ρ=true, verbose=false)
  x_approx = solve(S, b)
  @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
  @test x ≈ x_approx rtol=0.1

  for solver in ["splitBregman"]
    reg = Regularization("L1",1.e-3)
    solverInfo = SolverInfo(ComplexF64)
    S = createLinearSolver(solver,F; reg=reg,iterations=5,iterationsInner=40,
                                    ρ=1.0,solverInfo=solverInfo, normalizeReg=false)
    x_approx = solve(S, b)
    @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
    @test x ≈ x_approx rtol=0.1

    reg.λ *= length(b)/norm(b,1)
    S = createLinearSolver(solver,F; reg=reg,iterations=5,iterationsInner=40,
                            ρ=1.0,solverInfo=solverInfo, normalizeReg=true)
    x_approx = solve(S, b)
    @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
    @test x ≈ x_approx rtol=0.1
  end

  for solver in ["splitBregman"]
    reg = Regularization("L1",1.e-3)
    solverInfo = SolverInfo(ComplexF32)
    S = createLinearSolver(solver,F; reg=reg,iterations=5,iterationsInner=40,
                                    ρ=1.0,solverInfo=solverInfo, normalizeReg=false)
    x_approx = solve(S, b)
    @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
    @test x ≈ x_approx rtol=0.1

    reg.λ *= length(b)/norm(b,1)
    S = createLinearSolver(solver,F; reg=reg,iterations=5,iterationsInner=40,
                            ρ=1.0,solverInfo=solverInfo, normalizeReg=true)
    x_approx = solve(S, b)
    @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
    @test x ≈ x_approx rtol=0.1
  end

  for solver in ["primaldualsolver"]
    reg = [Regularization("L1",1.e-4), Regularization("TV",1.e-4)]
    solverInfo = SolverInfo(Float64)
    FR = [real.(F./norm(F)); imag.(F./norm(F))]
    bR = [real.(b./norm(F)); imag.(b./norm(F))]
    S = createLinearSolver(solver,FR; reg=reg, regName=["L1","TV"], iterations=1000, solverInfo=solverInfo)
    x_approx = solve(S, bR)
    @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
    @test x ≈ x_approx rtol=0.1
  end

  for solver in ["primaldualsolver"]
    reg = [Regularization("L1",1.e-4), Regularization("TV",1.e-4)]
    solverInfo = SolverInfo(Float32)
    FR = [real.(F./norm(F)); imag.(F./norm(F))]
    bR = [real.(b./norm(F)); imag.(b./norm(F))]
    S = createLinearSolver(solver,FR; reg=reg, regName=["L1","TV"], iterations=1000, solverInfo=solverInfo)
    x_approx = solve(S, bR)
    @info "Testing solver $solver ...: relative error = $(norm(x - x_approx) / norm(x))"
    @test x ≈ x_approx rtol=0.1
  end

end


text:
## -- Import and extend Base.display for various custom types

    import Base.display

    # Custom pretty printing for named tuples as datasets
    # TODO: possibly avoid blatant type piracy 🏴☠ 🏴☠ 🏴☠
    function display(x::NamedTuple)
        i = 1
        println("NamedTuple with $(length(keys(x))) elements:")
        l = max(length.(string.(keys(x)))...)
        for s in keys(x)
            t = typeof(x[s])
            sp = " "^(l-length(string(s)))
            print("  $s$sp  = $t")
            if t<:Number
                print("\t$(x[s])")
            elseif t<:AbstractRange
                print("\t$(x[s])")
            elseif t<:AbstractArray
                print(size(x[s]))
                if length(x[s]) < 2
                    print("\t[$(x[s])]")
                else
                    print("\t[$(x[s][1]) ... $(x[s][end])]")
                end
            elseif t<:AbstractString
                if length(x[s]) < 50
                    print("\t\"$(x[s])\"")
                else
                    print("\t\"$(x[s][1:50])...")
                end
            end
            print("\n")
            i += 1
            if i > 222
                print(".\n.\n.\n")
                break
            end
        end
    end


text:
# this is a work-in-progress test involving a picture of a duck

using Pixell, Plots
duck = read_map("data/duck.fits")
shape, wcs = size(duck), duck.wcs
plot(duck)

##
alm_duck = map2alm(duck)
nm = Pixell.alm2map(alm_duck, size(duck), duck.wcs)
plot(nm)

##
right_duck = duck[end:-1:begin, :]

plot(right_duck)

duckQU = cat(duck, right_duck; dims=3)
write_map("data/duckQU.fits", duckQU)

##
almDuck_QU = map2alm(QU)
newduck = alm2map(almDuck_QU, shape, wcs)
plot(newduck[1] .- duck)

##

pyQU = pycall(enmap.read_map, PyObject, "data/duckQU.fits")

myalm = pycall(curvedsky.map2alm, PyObject, pyQU, lmax=almDuck_QU[1].lmax, spin=2)
pyE = collect(myalm)[1,:][1]
pyB = collect(myalm)[2,:][1]

myalm2 = [pyE pyB]

##
almDuck_QU[1].alm
plot(abs.(almDuck_QU[1].alm), yscale=:log10, ylim=(1e-10,1e-1), lw=0.1)
plot!(abs.(pyE), yscale=:log10, ylim=(1e-10,1e-1), lw=0.1)

##
newpymap = pycall(enmap.enmap, PyObject, pyQU, pyQU.wcs)
# pyQUalm = collect(cat(almDuck_QU[1].alm,almDuck_QU[2].alm,dims=2)')
pyQUduck = pycall(curvedsky.alm2map, PyObject, PyReverseDims(myalm2), newpymap, spin=2)

##
l2 = py"lambda x : x[1,:,:]"
heatmap(l2(pyQUduck))

##
# heatmap(duck.data .- nm.data, clim=(-1e-5,1e-5))
#
# imap = pycall(enmap.read_map, PyObject, "test/data/duck.fits")
# pyshape = imap.shape
# pywcs = imap.wcs
# pya = curvedsky.map2alm(imap, lmax=alm_duck.lmax)
# newpymap = pycall(enmap.enmap, PyObject, imap,  imap.wcs)
# pycall(curvedsky.alm2map, PyObject, pya, newpymap, spin=0)
# heatmap(PyArray(imap)' - PyArray(newpymap)', clim=(-1e-5,1e-5))


text:
using StochasticDiffEq, LinearAlgebra, SparseArrays, Random, LinearSolve, Test
using StochasticDiffEq.OrdinaryDiffEq: WOperator, set_gamma!, calc_W!

@testset "Derivative Utilities" begin
  @testset "calc_W!" begin
    A = [-1.0 0.0; 0.0 -0.5]; σ = [0.9 0.0; 0.0 0.8]
    mm = [2.0 0.0; 0.0 1.0]
    u0 = [1.0, 1.0]; tmp = zeros(2)
    tspan = (0.0,1.0); dt = 0.01; dtgamma = 0.5dt
    concrete_W = -mm + dtgamma * A

    # Out-of-place
    _f = (u,p,t) -> A*u; _g = (u,p,t) -> σ*u
    fun = SDEFunction(_f, _g;
                      mass_matrix=mm,
                      jac=(u,p,t) -> A)
    prob = SDEProblem(fun, _g, u0, tspan)
    integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt)
    W = integrator.cache.nlsolver.cache.W
    calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, false)
    @test convert(AbstractMatrix, W) ≈ concrete_W
    @test W \ u0 ≈ concrete_W \ u0

    # In-place
    _f = (du,u,p,t) -> mul!(du,A,u); _g = (du,u,p,t) -> mul!(du,σ,u)
    fun = SDEFunction(_f, _g;
                      mass_matrix=mm,
                      jac_prototype=DiffEqArrayOperator(A))
    prob = SDEProblem(fun, _g, u0, tspan)
    integrator = init(prob, ImplicitEM(theta=1); adaptive=false, dt=dt)
    W = integrator.cache.nlsolver.cache.W
    calc_W!(W, integrator, integrator.cache.nlsolver, integrator.cache, dtgamma, false)

    # Did not update because it's an array operator
    # We don't want to build Jacobians when we have operators!
    @test convert(AbstractMatrix, integrator.cache.nlsolver.cache.W) != concrete_W
    ldiv!(tmp, lu!(integrator.cache.nlsolver.cache.W), u0); @test tmp != concrete_W \ u0

    # But jacobian2W! will update the cache
    StochasticDiffEq.OrdinaryDiffEq.jacobian2W!(integrator.cache.nlsolver.cache.W._concrete_form, mm, dtgamma, integrator.cache.nlsolver.cache.W.J.A, false)
    @test convert(AbstractMatrix, integrator.cache.nlsolver.cache.W) == concrete_W
    ldiv!(tmp, lu!(integrator.cache.nlsolver.cache.W), u0); @test tmp == concrete_W \ u0
  end

  @testset "Implicit solver with lazy W" begin
    A = sparse([-1.0 0.0; 0.0 -0.5]); σ = sparse([0.9 0.0; 0.0 0.8])
    mm = sparse([2.0 0.0; 0.0 1.0])
    u0 = [1.0, 1.0]; tspan = (0.0,1.0)

    _f = (u,p,t) -> t*(A*u); _f_ip = (du,u,p,t) -> lmul!(t,mul!(du,A,u))
    _g = (u,p,t) -> σ*u; _g_ip = (du,u,p,t) -> mul!(du,σ,u)
    prob1 = SDEProblem(SDEFunction(_f, _g; mass_matrix=mm), _g, u0, tspan)
    prob2 = SDEProblem(SDEFunction(_f, _g; mass_matrix=mm, jac=(u,p,t) -> t*A), _g, u0, tspan)
    prob1_ip = SDEProblem(SDEFunction(_f_ip, _g_ip; mass_matrix=mm), _g_ip, u0, tspan)
    jac_prototype=DiffEqArrayOperator(similar(A); update_func=(J,u,p,t) -> (J .= t .* A; J))
    prob2_ip = SDEProblem(SDEFunction(_f_ip, _g_ip; mass_matrix=mm, jac_prototype=jac_prototype), _g_ip, u0, tspan)

    for Alg in [ImplicitEM, ISSEM]
      println(Alg)
      Random.seed!(0); sol1 = solve(prob1, Alg(theta=1); adaptive=false, dt=0.01)
      Random.seed!(0); sol2 = solve(prob2, Alg(theta=1); adaptive=false, dt=0.01)
      @test sol1(1.0) ≈ sol2(1.0) rtol=1e-4
      Random.seed!(0); sol1_ip = solve(prob1_ip, Alg(theta=1); adaptive=false, dt=0.01)
      Random.seed!(0); sol2_ip = solve(prob2_ip, Alg(theta=1); adaptive=false, dt=0.01)
      @test sol1_ip(1.0) ≈ sol2_ip(1.0) rtol=1e-4
    end

    σ = 1.0
    _g = (u,p,t) -> σ; _g_ip = (du,u,p,t) -> (du .= σ)
    prob1 = SDEProblem(SDEFunction(_f, _g), _g, u0, tspan)
    prob2 = SDEProblem(SDEFunction(_f, _g; jac=(u,p,t) -> t*A), _g, u0, tspan)
    prob1_ip = SDEProblem(SDEFunction(_f_ip, _g_ip), _g_ip, u0, tspan)
    prob2_ip = SDEProblem(SDEFunction(_f_ip, _g_ip; jac_prototype=jac_prototype), _g_ip, u0, tspan)

    println(SKenCarp)
    Random.seed!(0); sol1 = solve(prob1, SKenCarp(); adaptive=false, dt=0.01)
    Random.seed!(0); sol2 = solve(prob2, SKenCarp(); adaptive=false, dt=0.01)
    @test sol1(1.0) ≈ sol2(1.0) rtol=1e-4
    Random.seed!(0); sol1_ip = solve(prob1_ip, SKenCarp(); adaptive=false, dt=0.01)
    Random.seed!(0); sol2_ip = solve(prob2_ip, SKenCarp(); adaptive=false, dt=0.01)
    @test sol1_ip(1.0) ≈ sol2_ip(1.0) rtol=1e-3
  end
end

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq

max_issues_repo_name: JuliaDiffEq/StochasticDiffEq

max_forks_repo_name: JuliaDiffEq/StochasticDiffEq