text:
import FITSIO
import JLD
import Optim

import ..Log
using ..Model
import ..SDSSIO
import ..Infer
import ..SDSSIO: RunCamcolField
import ..PSF

using ..DeterministicVI

function union_find!(i, components_tree)
    root = i
    while components_tree[i] != i
        i = components_tree[i]
    end

    while root != i
        root2 = components_tree[root]
        components_tree[root] = i
        root = root2
    end
    return i
end

"""
Computes connected components given a conflict graph, and the range of sources
to process. Writes to the components argument.

- source_start - start source
- end_source - end source (inclusively processed)
- sources - the actual source ids to find connected components for (indexed by source_start and source_end)
- neighbor_map - conflict graph of sources
- components_tree - the components array to do union find on.
- components - Dict{Int64, Vector{Int64}} dictionary from cc_id to target_source index in that component
- src_to_target_index - a mapping from light source id -> index within the target_sources array. This is required
                        since we need to write the index within target sources to the output, rather than the source id itself.

Returns:
- Nothing
"""
function compute_connected_components(source_start, source_end, sources, neighbor_map,
                                      components_tree, components, src_to_target_index)
    # Construct reverse map from source_id -> index in source. Only do this for
    # the specified portion we are computing connected components.
    source_to_index = Dict{Int64, Int64}()
    index_to_source = Dict{Int64, Int64}()
    for i = source_start:source_end
        source_to_index[sources[i]] = i-source_start+1
        index_to_source[i-source_start+1] = sources[i]
    end

    for i = source_start:source_end
        components_tree[i-source_start+1] = i-source_start+1
    end

    # Run union find algorithm to find connected components.
    for i = source_start:source_end
        target = union_find!(i-source_start+1, components_tree)
        for neighbor in neighbor_map[index_to_source[i-source_start+1]]
            if haskey(source_to_index, neighbor)
                neighbor_idx = source_to_index[neighbor]
                component_source = union_find!(neighbor_idx, components_tree)
                components_tree[component_source] = target
            end

        end
    end

    # Write to components dictionary.
    for i = source_start:source_end
        index = union_find!(i-source_start+1, components_tree)
        if !haskey(components, index)
            components[index] = Vector{Int64}()
        end
        push!(components[index], src_to_target_index[index_to_source[i-source_start+1]])
    end
end

"""
Partitions sources via the cyclades algorithm.

- nprocthreads - number of threads to which to distribute sources
- target_sources - array of target sources. Elements should match keys of neighbor_map.
- neighbor_map - graph of connections of sources

Returns:
- An array of vectors representing the workload of each thread ([thread][batch][sources(indices)])
"""
function partition_cyclades(nprocthreads, target_sources, neighbor_map; batch_size=60)
    Log.info("Starting Cyclades partitioning...")
    tic()

    n_sources = length(target_sources)
    n_total_batches = convert(Int64, ceil(n_sources / batch_size))

    # Construct src_to_indx map
    src_to_indx = Dict{Int64, Int64}()
    for i=1:n_sources
        src_to_indx[target_sources[i]] = i
    end

    # The final workload distribution.
    thread_sources_assignment = Array(Vector{Vector{Int64}}, nprocthreads)
    for thread = 1:nprocthreads
        thread_sources_assignment[thread] = Vector{Vector{Int64}}(n_total_batches)
        for batch = 1:n_total_batches
            thread_sources_assignment[thread][batch] = Vector{Int64}()
        end
    end

    # First shuffle the sources. Note Cyclades is serially equivalent
    # to this permutation of sources.
    sources = [x for x in keys(neighbor_map)]
    shuffle!(sources)

    # Process batch_size number of sources at a time.
    # TODO max - parallelize everything below (particularly CC computation)

    # The components tree is for union-find.
    components_tree = [Array(Int64, batch_size) for i=1:nprocthreads]

    # We have n_total_batches components, where each component is a dictionary
    # from the component id, to a list of sources in that component
    components = [Dict{Int64, Vector{Int64}}() for i=1:n_total_batches]
    sources_in_components = 0
    for (cur_batch, source_index) in enumerate(collect(1:batch_size:n_sources))
        source_start = source_index
        source_end = min(source_index+batch_size-1, n_sources)
        @assert source_end - source_start + 1 <= batch_size

        # TODO max - parallelize this
        compute_connected_components(source_start, source_end, sources, neighbor_map,
                                     components_tree[1],
                                     components[cur_batch],
                                     src_to_indx)
    end

    assigned_sources = 0
    # Load balance the connected components within each batch into thread_sources_assignment.
    for (cur_batch, cur_batch_component) in enumerate(components)
        # Priority queue for load balancing.
        pqueue = Base.Collections.PriorityQueue([i for i=1:nprocthreads],
                                                [0 for i=1:nprocthreads],
                                                Base.Order.Forward)

        # Assign non-conflicting group of sources to different threads
        for (component_group_id, sources_of_component) in cur_batch_component
            least_loaded_thread = Base.Collections.peek(pqueue)[1]
            for source_to_assign in sources_of_component
                push!(thread_sources_assignment[least_loaded_thread][cur_batch], source_to_assign)
                assigned_sources += 1
            end
            pqueue[least_loaded_thread] += length(sources_of_component)
        end
    end

    Log.info("Cyclades - Assigned sources: $(assigned_sources) vs correct number of sources: $(n_sources)")
    @assert assigned_sources == n_sources
    Log.info("Cyclades - Number of batches: $(n_total_batches)")
    Log.info("Finished Cyclades partitioning.  Elapsed time: $(toq()) seconds")
    thread_sources_assignment
end

"""
Partitions sources across threads equally.

- nprocthreads - number of threads to which to distribute processing the sources.
- n_sources - the number of total sources to process.

Returns:
- An array of vectors representing the workload of each thread ([thread][batch][sources]).
  Note for partition equally, there is only 1 batch.

"""
function partition_equally(nprocthreads, n_sources)
    Log.info("Starting basic source partitioning...")
    tic()
    n_sources_per_thread = floor(Int64, n_sources / nprocthreads)
    thread_sources_assignment = Array(Vector{Vector{Int64}}, nprocthreads)
    n_sources_assigned = 0
    for thread = 1:nprocthreads
        start_source = (thread-1) * n_sources_per_thread
        end_source = thread * n_sources_per_thread
        if thread == nprocthreads
            end_source = n_sources
        end
        # Only 1 batch for partitioning equally
        thread_sources_assignment[thread] = Vector{Vector{Int64}}(1)
        thread_sources_assignment[thread][1] = Vector{Int64}(end_source-start_source)
        for source_assignment = start_source:end_source-1
            thread_sources_assignment[thread][1][source_assignment-start_source+1] = source_assignment+1
            n_sources_assigned += 1
        end
    end
    @assert n_sources_assigned == n_sources
    Log.info("Finished basic source partitioning. Elapsed time: $(toq()) seconds")
    thread_sources_assignment
end

"""
Like one_node_infer, uses multiple threads on one node to fit the Celeste
model over numerous iterations.

catalog - the catalog of light sources
target_sources - light sources to optimize
neighbor_map - ligh_source index -> neighbor light_source id

cyclades_partition - use the cyclades algorithm to partition into non conflicting batches for updates.
batch_size - size of a single batch of sources for updates
within_batch_shuffling - whether or not to process sources within a batch randomly
joint_inference_terminate - whether to terminate once sources seem to be stable
termination_percent - stop optimization once a certain percentage of sources have been optimized.

Returns:

- Vector of OptimizedSource results
"""
function one_node_joint_infer(catalog, target_sources, neighbor_map, images;
                              cyclades_partition=true,
                              batch_size=60,
                              within_batch_shuffling=true,
                              termination_percent=.95,
                              n_iters=10)
    # Seed random number generator to ensure the same results per run.
    srand(42)

    nprocthreads = nthreads()

    # Partition the sources
    n_sources = length(target_sources)
    Log.info("Optimizing $(n_sources) sources")
    if cyclades_partition
        # Convert neighbormap to cyclades map (map from source_id -> [neighbor source_ids])
        cyclades_neighbor_map = Dict{Int64, Vector{Int64}}()
        for (index, neighbors) in enumerate(neighbor_map)
            source_id = target_sources[index]
            cyclades_neighbor_map[source_id] = neighbors
        end
        thread_sources_assignment = partition_cyclades(nprocthreads, target_sources, cyclades_neighbor_map, batch_size=batch_size)
    else
        thread_sources_assignment = partition_equally(nprocthreads, n_sources)
    end

    Log.info("Done assigning sources to threads for processing")

    # Pre allocate elbo args variational params
    target_source_variational_params = Dict{Int64, Array{Float64}}()
    for target_source in target_sources
        cat = catalog[target_source]
        target_source_variational_params[target_source] = init_source(cat)
    end

    # Pre-allocate dictionary of elboargs, call it ea_vec.
    ea_vec = Array{ElboArgs}(n_sources)
    function initialize_elboargs_sources(sources)
#        nputs(dt_nodeid, "Thread $(Threads.threadid()) allocating mem for $(length(sources)) sources")
        for cur_source_index in sources
            entry_id = target_sources[cur_source_index]
            entry = catalog[target_sources[cur_source_index]]
            neighbor_ids = neighbor_map[cur_source_index]
            neighbors = catalog[neighbor_map[cur_source_index]]

            # TODO max: refactor this portion? It's reused in infer_source.
#            nputs(dt_nodeid, "Thread $(Threads.threadid()) allocating mem for source $(target_sources[cur_source_index]): objid=$(entry.objid)")
            cat_local = vcat([entry], neighbors)
            ids_local = vcat([entry_id], neighbor_ids)

            #vp = Vector{Float64}[init_source(ce) for ce in cat_local]
            vp = Vector{Float64}[haskey(target_source_variational_params, x) ?
                        target_source_variational_params[x] :
                        init_source(catalog[x]) for x in ids_local]
            patches = Infer.get_sky_patches(images, cat_local)
            ea = ElboArgs(images, vp, patches, [1])
            Infer.load_active_pixels!(ea.images, ea.patches)
            ea_vec[cur_source_index] = ea
        end
    end

    # Initialize elboargs in parallel
    tic()
    thread_initialize_sources_assignment = partition_equally(nprocthreads, n_sources)

    Threads.@threads for i=1:nprocthreads
        for batch = 1:length(thread_initialize_sources_assignment[i])
            initialize_elboargs_sources(thread_initialize_sources_assignment[i][batch])
        end
    end

    Log.info("Done preallocating array of elboargs. Elapsed time: $(toq())")

    # Keep track of which sources have converged.
    sources_converged = Dict{Int64, Bool}()
    for source in target_sources
        sources_converged[source] = false
    end
    n_sources_converged = 0
    n_sources_converged_lock = SpinLock()

    function should_optimize_source(src_indx)
        src_has_converged = sources_converged[target_sources[src_indx]]
        neighbors_have_converged = true
        for neighbor in neighbor_map[src_indx]
            if haskey(sources_converged, neighbor)
                neighbors_have_converged = neighbors_have_converged && sources_converged[neighbor]
            end
        end
        return !src_has_converged || !neighbors_have_converged
    end

    # Process partition of sources. Multiple threads call this function in parallel.
    function process_sources(source_assignment::Vector{Int64}, iter)
        try

            # Use a constant number of newton steps
            n_newton_steps = 5

            # Shuffle the source assignments within each batch of each process.
            # This is disabled by default because it ruins the deterministic outcome
            # required by the test cases.
            if within_batch_shuffling
                shuffle!(source_assignment)
            end
            for cur_source_indx in source_assignment
                # Optimize only if source has not converged or at least
                # one of its neighbors has not converged
                if should_optimize_source(cur_source_indx)
                    cur_entry = catalog[target_sources[cur_source_indx]]
                    iter_count, obj_value, max_x, r = DeterministicVI.maximize_f(
                        DeterministicVI.elbo,
                        ea_vec[cur_source_indx],
                        max_iters=n_newton_steps,
                        use_default_optim_params=true)
                    sources_converged[target_sources[cur_source_indx]] = Optim.converged(r)
                end

                # Maintain count of sources that have converged
                if sources_converged[target_sources[cur_source_indx]]
                    lock(n_sources_converged_lock)
                    n_sources_converged += 1
                    unlock(n_sources_converged_lock)
                end
            end
        catch ex
            if is_production_run || nthreads() > 1
                Log.error(string(ex))
            else
                rethrow(ex)
            end
        end
    end

    # Process sources in parallel using nprocthreads.
    tic()
    n_batches = length(thread_sources_assignment[1])
    for iter = 1:n_iters
        # Reset number of sources converged
        n_sources_converged = 0

        # Process every batch of every iteration. We do the batches on the outside
        # Since there is an implicit barrier after the inner threaded for loop below.
        # We want this barrier because there may be conflict _between_ Cyclades batches.
        for batch = 1:n_batches
            # Process every batch of every iteration with nprocthreads
            Threads.@threads for i = 1:nprocthreads
                process_sources(thread_sources_assignment[i][batch], iter)
            end
        end

        if n_sources_converged >= termination_percent * n_sources
            break
        end
    end
    Log.info("$(n_sources_converged) / $(n_sources) converged")
    Log.info("Done fitting elboargs. Elapsed time: $(toq())")

    # Return add results to vector
    results = OptimizedSource[]

    for i = 1:n_sources
        entry = catalog[target_sources[i]]
        result = OptimizedSource(entry.thing_id,
                                 entry.objid,
                                 entry.pos[1],
                                 entry.pos[2],
                                 ea_vec[i].vp[1])
        push!(results, result)
    end

    results
end


text:
# Example for basic 2D modeling:
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#

using SetIntersectionProjection
using JUDI, SeisIO, JOLI, PyPlot, ProximalOperators, IterativeSolvers, SlimOptim, MAT, Dierckx, SpecialFunctions, LinearAlgebra, SparseArrays, FFTW
using Statistics, Random

## setup for bas
mutable struct compgrid
  d :: Tuple
  n :: Tuple
end

function spg_tv_norm(d_sim,objective_function,q_true, q_init, model; maxiter =5,tv_norm = 0.5)

  # with constraints:
  options=PARSDMM_options()
  options.FL=Float32
  options=default_PARSDMM_options(options,options.FL)
  options.adjust_gamma           = true
  options.adjust_rho             = true
  options.adjust_feasibility_rho = true
  options.Blas_active            = true
  options.maxit                  = 1000
  options.feas_tol= 0.001
  options.obj_tol=0.001
  options.evol_rel_tol = 0.00001

  options.rho_ini=[1.0f0]

  set_zero_subnormals(true)
  BLAS.set_num_threads(2)
  FFTW.set_num_threads(2)
  options.parallel=false
  options.feasibility_only = false
  options.zero_ini_guess=true

  #select working precision
  if options.FL==Float64
    TF = Float64
    TI = Int64
  elseif options.FL==Float32
    TF = Float32
    TI = Int32
  end

  #----# make projector
  #comp_grid=compgrid((44.75,44.75),(model.n[1],model.n[2]))
  comp_grid=compgrid((1,1),(model.n[1],model.n[2]))

  constraint = Vector{SetIntersectionProjection.set_definitions}()

  #TV
  (TV,dummy1,dummy2,dummy3)=get_TD_operator(comp_grid,"TV",options.FL)
  m_min     = 0.0
  m_max     = tv_norm*norm(TV*vec(q_true),1)
  set_type  = "l1"
  TD_OP     = "TV"
  app_mode  = ("matrix","")
  custom_TD_OP = ([],false)
  push!(constraint, set_definitions(set_type,TD_OP,m_min,m_max,app_mode,custom_TD_OP))

  #print(0.5f0*norm(TV*vec(weights),1))

  #NEW
  (P_sub,TD_OP,set_Prop) = setup_constraints(constraint,comp_grid,options.FL)
  (TD_OP,AtA,l,y)        = PARSDMM_precompute_distribute(TD_OP,set_Prop,comp_grid,options)
  options.rho_ini        = ones(length(TD_OP))*10.0

  proj_intersection = x_new-> PARSDMM(x_new,AtA,TD_OP,set_Prop,P_sub,comp_grid,options)

  function prj!(input)
      (x, dummy1, dummy2, dymmy3) = proj_intersection(input)
      return x
  end

  #proj_q_adj = prj!(vec(q_norm))
  #imshow(reshape(proj_q_adj,n))

  # Source localization with SPG
  options_spg            = spg_options(verbose=3, maxIter=maxiter, memory=5)
  options_spg.progTol	   = 1e-60

  #options_spg = spg_options(verbose=3, maxIter=maxiter, memory=5, suffDec=1f-6)
  #options_spg.progTol=1e-60
  #options_spg.testOpt=false
  #options_spg.interp=2

  #ProjBound(x)       = x

  #print(length(w_init))
  q_spg, fsave, funEvals = minConf_SPG(objective_function, vec(q_init), prj!, options_spg)

  return q_spg
end

ext: jl


text:
# This example highlights the use of a linear minimization oracle
# using an LP solver defined in MathOptInterface
# we compare the performance of the two LMOs, in- and out of place
#
# to get accurate timings it is important to run twice so that the
# compile time of Julia for the first run is not tainting the results

include("activate.jl")

using JSON
using LaTeXStrings
results = JSON.Parser.parsefile("lcg_expensive_data.json")

ref_optimum = results["reference_BCG_primal"]

iteration_list = [
    [x[1] + 1 for x in results["FW"]],
    [x[1] + 1 for x in results["LCG"]],
    [x[1] + 1 for x in results["BLCG"]],
    [x[1] + 1 for x in results["LAFW"]],
    [x[1] + 1 for x in results["BCG"]],
]
time_list = [
    [x[5] for x in results["FW"]],
    [x[5] for x in results["LCG"]],
    [x[5] for x in results["BLCG"]],
    [x[5] for x in results["LAFW"]],
    [x[5] for x in results["BCG"]],
]
primal_gap_list = [
    [x[2] - ref_optimum for x in results["FW"]],
    [x[2] - ref_optimum for x in results["LCG"]],
    [x[2] - ref_optimum for x in results["BLCG"]],
    [x[2] - ref_optimum for x in results["LAFW"]],
    [x[2] - ref_optimum for x in results["BCG"]],
]
dual_gap_list = [
    [x[4] for x in results["FW"]],
    [x[4] for x in results["LCG"]],
    [x[4] for x in results["BLCG"]],
    [x[4] for x in results["LAFW"]],
    [x[4] for x in results["BCG"]],
]
label = [L"\textrm{FW}", L"\textrm{L-CG}", L"\textrm{BL-CG}", L"\textrm{L-AFW}", L"\textrm{BCG}"]

FrankWolfe.plot_results(
    [primal_gap_list, primal_gap_list, dual_gap_list, dual_gap_list],
    [iteration_list, time_list, iteration_list, time_list],
    label,
    [L"\textrm{Iteration}", L"\textrm{Time}", L"\textrm{Iteration}", L"\textrm{Time}"],
    [L"\textrm{Primal Gap}", L"\textrm{Primal Gap}", L"\textrm{Dual Gap}", L"\textrm{Dual Gap}"],
    xscalelog=[:log, :identity, :log, :identity],
    legend_position=[:bottomleft, nothing, nothing, nothing],
    filename="lcg_expensive.pdf",
)


text:
using Printf

function output_result(stepnum, Qbase, cellxmax, cellymax, cellzmax, out_file_front, out_ext, out_dir, Rd)

    stepnum = string(stepnum)
    while length(stepnum) < 6
        stepnum = "0"*stepnum
    end

    fff = out_dir*"/" * out_file_front*stepnum*out_ext
    open(fff,"w") do f
        write(f,"result:rho[kg/m^3], u[m/s], v[m/s], w[m/s], p[Pa], T[K]\n")
        for i in 2:cellxmax-1
            for j in 2:cellymax-1
                for k in 2:cellzmax-1
                    a1 = @sprintf("%8.8e", Qbase[i,j,k,1])
                    a2 = @sprintf("%8.8e", Qbase[i,j,k,2])
                    a3 = @sprintf("%8.8e", Qbase[i,j,k,3])
                    a4 = @sprintf("%8.8e", Qbase[i,j,k,4])
                    a5 = @sprintf("%8.8e", Qbase[i,j,k,5])
                    T = Qbase[i,j,k,5]/(Qbase[i,j,k,1]*Rd)
                    a6 = @sprintf("%8.8e", T)
                    write(f, a1*" "*a2*" "*a3*" "*a4*" "*a5*" "*a6*"\n")
                end
            end
        end
    end
    println("\nwrite "*fff)
end

function reset_write(fwrite)
    open( fwrite, "w" ) do f
    end
end

function output_physicaltime(fwrite, t, dt)
    a = string(t)
    b = @sprintf("%8.8f", t*dt)
    open(fwrite, "a") do f
        write(f, "\n")
        write(f, "step = "*a)
        write(f, "     ")
        write(f, "time[s] = "*b)
        write(f, "\n")
        write(f, "\n")
    end
end

function output_innertime(fwrite, tau, norm2)
    a = string(tau)
    open(fwrite, "a") do f
        write(f, a)
        write(f, " ")
        for l in 1:5
            a = @sprintf("%8.8e", norm2[l])
            write(f, a)
            write(f, " ")
        end
        write(f, "\n")
    end
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
#include("quantumoptics.jl")
#include("ode.jl")

using quantumoptics
#using ode

function profile_basis()
	basis0 = GenericBasis([4,5])
	basis1 = GenericBasis([3,7])
	basis = bases.compose(basis0, basis1)
	#basis = basis0

	println(bases.multiplicable(basis, basis))
	function f(N)
		for i=1:N
			#basis==basis
			bases.multiplicable(basis, basis)
		end
	end
	f(1)
	@time f(10000)
end

function profile_mul()
	basis = FockBasis(3)
	op2 = destroy(basis)
	op1 = create(basis)
	tmp = Operator(basis,basis)
	tmp = op1*op2
	@time tmp = op1*op2
	println(tmp)

	inplacearithmetic.mul!(op1, op2, tmp)
	println(tmp)
	@time inplacearithmetic.mul!(op1, op2, tmp)
	inplacearithmetic.mul!(op1.data, op2.data, tmp.data)
	@time inplacearithmetic.mul!(op1.data, op2.data, tmp.data)
end

function profile_add()
	basis = FockBasis(400)
	op1 = destroy(basis)
	op2 = create(basis)
	function f(N)
		for i=1:N
			op1 += op2
		end
	end
	f(1)
	@time f(100)

	op1 = destroy(basis)
	op2 = create(basis)
	result = Operator(basis)
	function f2(N)
		for i=1:N
			inplacearithmetic.add!(op1,op2,result)
		end
	end
	f2(1)
	@time f2(100)
end

function profile_dmaster()
	basis_cavity = FockBasis(2^11)
	a = destroy(basis_cavity)
 	x0 = basis_ket(basis_cavity, 3) + basis_ket(basis_cavity, 5)
 	rho0 = tensor(x0, dagger(x0))
	H = 1*(a + dagger(a))
	J = [a]
    Jdagger = [dagger(j) for j=J]
    Heff = -1im*H -0.5*dagger(a)*a
    Heff_dagger = dagger(Heff)
    tmp1 = Operator(rho0.basis_l, rho0.basis_r)
    tmp2 = Operator(rho0.basis_l, rho0.basis_r)
    drho = Operator(rho0.basis_l, rho0.basis_l)

    rho0_data = rho0.data
    Heff_data = Heff.data
    J_data = [j.data for j=J]
    Jdagger_data = [j.data for j=Jdagger]
    drho_data = drho.data
    tmp_data = tmp1.data
    function f(N)
        for i=1:N
            quantumoptics.timeevolution_simple.dmaster(rho0, H, J, Jdagger)
        end
    end
    function f2(N)
        for i=1:N
            quantumoptics.timeevolution.dmaster(rho0, H, J, Jdagger, [drho, tmp1, tmp2])
        end
    end
    function f3(N)
        for i=1:N
            quantumoptics.timeevolution.dmaster_Heff(rho0_data, Heff_data, Heff_dagger.data, J_data, Jdagger_data, drho_data, tmp_data)
        end
    end
    function f4(N)
        for i=1:N
            quantumoptics.timeevolution.dmaster_H(rho0_data, H.data, J_data, Jdagger_data, drho_data, tmp_data)
        end
    end
    function f5(N)
        for i=1:N
            quantumoptics.timeevolution.dmaster_nh(rho0_data, Heff.data, Heff_dagger.data, J, Jdagger, drho_data, tmp_data)
        end
    end
    drho = quantumoptics.timeevolution_simple.dmaster(rho0, H, J, Jdagger)
    quantumoptics.timeevolution.dmaster_H(rho0_data, H.data, J_data, Jdagger_data, drho_data, tmp_data)
    drho1 = drho_data*1
    quantumoptics.timeevolution.dmaster_Heff(rho0_data, Heff, Heff_dagger, J, Jdagger, drho_data, tmp_data)
    drho2 = drho_data*1
    println(norm(drho2-drho.data))
    println(norm(drho1-drho.data))
    println(norm(drho1-drho2))
    #error()
    @time 1

    f(1)
    @time f(100)
    f2(1)
    @time f2(100)
    f3(1)
    @time f3(100)
    f4(1)
    @time f4(100)
    f5(1)
    @time f5(100)
    #@profile f4(10000)
end

function profile_master()
	basis_cavity = FockBasis(2^6)
	a = destroy(basis_cavity)
	x0 = basis_ket(basis_cavity, 3)
	rho0 = tensor(x0, dagger(x0))
	H = 0.5*(a + dagger(a)) + a*a + dagger(a)*dagger(a)
    #H = H + 0.1*Operator(basis_cavity, rand(2^5,2^5))
    H = H+H'
    J = [a]
	#J = [0.01*(a+0.01*Operator(basis_cavity, rand(2^5,2^5)))]
	Jdagger = [dagger(J[1])]
	T = [0.,1.]

    H_eff = H -1im*0.5*Jdagger[1]*J[1]
    H_eff_dagger = dagger(H_eff)
    tmpop = Operator(rho0.basis_l, rho0.basis_r)
    tmp = tmpop.data
    J_data = [j.data for j=J]
    Jdagger_data = [j.data for j=Jdagger]
    tmps = Matrix{Complex128}[]
    for i=length(tmps):11
        push!(tmps, Operator(rho0.basis_l, rho0.basis_r).data)
    end
    tmps = [H_eff.data*1 for i=1:11]
    tmps_op = Operator[]
    for i=1:14
        push!(tmps_op, Operator(rho0.basis_l, rho0.basis_r))
    end

    @profile 1
    @time 1

	tout, xout = quantumoptics.timeevolution_simple.master(T, rho0, H, J)
	@time quantumoptics.timeevolution_simple.master(T, rho0, H, J)

	t, x0 = quantumoptics.timeevolution.master(T, rho0, H, J, Jdagger=Jdagger, tmp=tmpop, tmps_ode=tmps)
    @time quantumoptics.timeevolution.master(T, rho0, H, J, Jdagger=Jdagger, tmp=tmpop, tmps_ode=tmps)

    t, x1 = quantumoptics.timeevolution.master_nh(T, rho0, H_eff, J, Hdagger=H_eff_dagger, Jdagger=Jdagger, tmp=tmpop, tmps_ode=tmps)
    @time quantumoptics.timeevolution.master_nh(T, rho0, H_eff, J, Hdagger=H_eff_dagger, Jdagger=Jdagger, tmp=tmpop, tmps_ode=tmps)
    #@profile quantumoptics.timeevolution.master_nh(T, rho0, H_eff, J, Hdagger=H_eff_dagger, Jdagger=Jdagger, tmp=tmpop, tmps_ode=tmps)
    println(norm(x1[end]-x0[end], Inf))
	println(norm(x1[end]-xout[end], Inf))
    println(norm(x0[end]-xout[end], Inf))
	#println(xout2[end])
end

Profile.clear()
#profile_basis()
#profile_mul()
#profile_add()
#println(methods(operators.zero!))
#profile_dmaster()
profile_master()
# using ProfileView
# ProfileView.view()
0


text:
function exponential(x::Vector)
    return exp((2.0 - x[1])^2) + exp((3.0 - x[2])^2)
end

function exponential_gradient!(storage::Vector, x::Vector)
    storage[1] = -2.0 * (2.0 - x[1]) * exp((2.0 - x[1])^2)
    storage[2] = -2.0 * (3.0 - x[2]) * exp((3.0 - x[2])^2)
end

function exponential_hessian!(storage::Matrix, x::Vector)
    storage[1, 1] = 2.0 * exp((2.0 - x[1])^2) * (2.0 * x[1]^2 - 8.0 * x[1] + 9)
    storage[1, 2] = 0.0
    storage[2, 1] = 0.0
    storage[2, 2] = 2.0 * exp((3.0 - x[2])^2) * (2.0 * x[2]^2 - 12.0 * x[2] + 19)
end

function nlsolvers_exponential(x, gx=nothing, hx=nothing)
	fx = exponential(x)
	if gx !== nothing
		exponential_gradient!(gx, x)
	end
	if hx !== nothing
		exponential_hessian!(hx, x)
	end
	objective_return(fx, gx, hx)
end
x0 =[0.0, 0.0]
xopt = [2.0, 3.0]
fxopt = exponential([2.0, 3.0])

##########################################################################
###
### Fletcher-Powell
###
### From [2]
### Source: A rapidly convergent descent method for minimization
###         Fletcher & Powell
##########################################################################

function fletcher_powell(x::Vector)
    function theta(x::Vector)
        if x[1] > 0
            return atan(x[2] / x[1]) / (2.0 * pi)
        else
            return (pi + atan(x[2] / x[1])) / (2.0 * pi)
        end
    end

    return 100.0 * ((x[3] - 10.0 * theta(x))^2 +
        (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2
end

function fletcher_powell_gradient!(storage::Vector, x::Vector)
    function theta(x::Vector)
        if x[1] > 0
            return atan(x[2] / x[1]) / (2.0 * pi)
        else
            return (pi + atan(x[2] / x[1])) / (2.0 * pi)
        end
    end

    if ( x[1]^2 + x[2]^2 == 0 )
        dtdx1 = 0;
        dtdx2 = 0;
    else
        dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
        dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
    end

    storage[1] = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
        200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 );
    storage[2] = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
        200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 );
    storage[3] =  200.0*(x[3]-10.0*theta(x)) + 2.0*x[3];
end

function fletcher_powell_fun_gradient!(storage::Vector, x::Vector)
    function theta(x::Vector)
        if x[1] > 0
            return atan(x[2] / x[1]) / (2.0 * pi)
        else
            return (pi + atan(x[2] / x[1])) / (2.0 * pi)
        end
    end

    if ( x[1]^2 + x[2]^2 == zero(eltype(x)) )
        dtdx1 = zero(eltype(x))
        dtdx2 = zero(eltype(x))
    else
        dtdx1 = - x[2] / ( 2.0 * pi * ( x[1]^2 + x[2]^2 ) )
        dtdx2 =   x[1] / ( 2.0 * pi * ( x[1]^2 + x[2]^2 ) )
    end

    storage[1] = -2000.0*(x[3]-10.0*theta(x))*dtdx1 +
        200.0*(sqrt(x[1]^2+x[2]^2)-1.0)*x[1]/sqrt( x[1]^2+x[2]^2 )
    storage[2] = -2000.0*(x[3]-10.0*theta(x))*dtdx2 +
        200.0*(sqrt(x[1]^2+x[2]^2)-1.0)*x[2]/sqrt( x[1]^2+x[2]^2 )
    storage[3] =  200.0*(x[3]-10.0*theta(x)) + 2.0*x[3]
    return 100.0 * ((x[3] - 10.0 * theta(x))^2 +
        (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2
end

function nlsolvers_fletcher_powell(x, gx=nothing, hx=nothing)
	fx = fletcher_powell(x)
	if gx !== nothing
		fletcher_powell_gradient!(gx, x)
	end
	if hx !== nothing
		fletcher_powell_hessian!(hx, x)
	end
	objective_return(fx, gx, hx)
end

x0=[-1.0, 0.0, 0.0] # Same as in source
xopt=[1.0, 0.0, 0.0]
0.0


text:
using PyPlot, ElasticWave
nz = 133; nx = 203; ext= 20; iflag = 1;
dx = 5. ; dz = 5. ; dt = 5e-4; tmax = 0.5; f0=30.;
vp =4000.*ones(nz,nx);  vs=2000.*ones(nz,nx); rho=2.5*ones(nz,nx)*1e-4;
fidMtx = CreateFidMtx(nz, nx, ext, iflag, vp, vs, rho, dz, dx, dt, f0);

if iflag == 1
   Nz = nz + 2*ext
elseif iflag == 2
   Nz = nz +   ext
end
Nx = nx + 2*ext

root = homedir();
path_pv = join([root "/Desktop/pv.bin"])
isx = collect(1:2:nx); isz = ones(Int64,length(isx)); ot = zeros(length(isx));
flags = vec([false false true true false]);
srcs  = InitMultiSources(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);
MultiStepForward(path_pv, srcs, fidMtx, dz, dx)
(vxx, vxz, vzx, vzz) = ReadPv(path_pv, 300); p = vxx+vzz; s=vxz-vzx;
SeisPlot(reshape(p, Nz, Nx)); SeisPlot(reshape(s, Nz, Nx));

path_wfd = join([root "/Desktop/wfd.bin"])
MultiStepForward(path_wfd, srcs, fidMtx)
wfd = ReadWfd(path_wfd, 300);
SeisPlot(reshape(wfd.Vx, Nz, Nx))

dm = zeros(Nz, Nx); du = zeros(Nz, Nx); dm[floor(Int64, Nz/2), 1:end] = 500.; du[floor(Int64, Nz/2), 1:end] = 500.;
dm = vec(dm); du = vec(du);
irx = collect(1:2:nx); irz = 2*ones(Int64, length(irx));
shotv = MultiStepForward(irz, irx, path_pv,  dm, du, fidMtx)
SeisPlot(shotv.Vz); SeisPlot(shotv.Vx);

(dm1, du1) = MultiStepAdjoint(shotv, fidMtx, path_pv);
SeisPlot(reshape(dm1, Nz, Nx)); SeisPlot(reshape(du1, Nz, Nx));

# ==============================================================================
using PyPlot, ElasticWave

function main()
    nz = 177; nx = 203; ext= 30;   iflag = 1;
    dx = 5. ; dz = 5. ; dt = 4e-4; tmax = 0.5; f0=30.;
    vp =3500.*ones(nz,nx);  vs=1750.*ones(nz,nx); rho=2.5*ones(nz,nx);
    vp[77:end,:] = 5000.; vs[77:end,:] = 2500.
    fidMtx = CreateFidMtx(nz, nx, ext, iflag, vp, vs, rho, dz, dx, dt, f0);

    isx = collect(1:5:nx); isz = ones(Int64,length(isx)); ot = zeros(length(isx));
    flags = vec([false false true true false]);
    srcs  = InitMultiSources(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);

    irx   = collect(1:1:nx); irz = 1*ones(Int64, length(irx));
    shotv = MultiStepForward(irz, irx, srcs, fidMtx)
    imshow(shotv.Vx, aspect=0.2, vmax=maximum(abs(shotv.Vx)), vmin=-maximum(abs(shotv.Vx)), cmap="PuOr"); figure();
    imshow(shotv.Vz, aspect=0.2, vmax=maximum(abs(shotv.Vz)), vmin=-maximum(abs(shotv.Vz)), cmap="PuOr")
    root = homedir()
    path_shotv = join([root "/Desktop/shotv.bin"])
    WriteShotV(path_shotv, shotv);

    shotv_FA = removeFirstArrival();
    shotv.Vx = shotv.Vx - shotv_FA.Vx;
    shotv.Vz = shotv.Vz - shotv_FA.Vz;
    figure(); imshow(shotv.Vx, aspect=0.2, vmax=maximum(abs(shotv.Vx)), vmin=-maximum(abs(shotv.Vx)), cmap="PuOr");
    figure(); imshow(shotv.Vz, aspect=0.2, vmax=maximum(abs(shotv.Vz)), vmin=-maximum(abs(shotv.Vz)), cmap="PuOr");

    path_shotv_DFA = join([root "/Desktop/shotv_DFA.bin"])
    WriteShotV(path_shotv_DFA, shotv);

    path_pv = join([root "/Desktop/pv.bin"])
    MultiStepForward(path_pv, srcs, fidMtx, dz, dx)
    (dm, du) = MultiStepAdjoint(shotv, fidMtx, path_pv);
    figure(); a = quantile(abs(dm[:]),0.98); imshow(dm, vmax=a, vmin=-a, cmap="PuOr");
    figure(); a = quantile(abs(du[:]),0.98); imshow(du, vmax=a, vmin=-a, cmap="PuOr");

    vps = modSmooth(vp, 20); vss = modSmooth(vs,20);
    fidMtx_SM = CreateFidMtx(nz, nx, ext, iflag, vps, vss, rho, dz, dx, dt, f0);

    path_pvsm = join([root "/Desktop/pv_sm.bin"])
    MultiStepForward(path_pvsm, srcs, fidMtx_SM, dz, dx)
    (dm2, du2) = MultiStepAdjoint(shotv, fidMtx_SM, path_pvsm);
    figure(); a = quantile(abs(dm1[:]),0.98); imshow(dm1, vmax=a, vmin=-a, cmap="PuOr");
    figure(); a = quantile(abs(du1[:]),0.98); imshow(du1, vmax=a, vmin=-a, cmap="PuOr");

    figure(); a = quantile(abs(dm2[:]),0.98); imshow(dm2, vmax=a, vmin=-a, cmap="PuOr");
    figure(); a = quantile(abs(du2[:]),0.98); imshow(du2, vmax=a, vmin=-a, cmap="PuOr");
end

function removeFirstArrival()
    nz = 177; nx = 203; ext= 30;   iflag = 1;
    dx = 10. ; dz = 10. ; dt = 1e-3; tmax = 1.0; f0=30.;
    vp =3500.*ones(nz,nx);  vs=1750.*ones(nz,nx); rho=2.5*ones(nz,nx);
    fidMtx = CreateFidMtx(nz, nx, ext, iflag, vp, vs, rho, dz, dx, dt, f0);

    isx = collect(1:5:nx); isz = ones(Int64,length(isx)); ot = zeros(length(isx));
    flags = vec([false false true true false]);
    srcs  = InitMultiSources(isz, isx, nz, nx, ext, iflag, f0, ot, dt, flags);

    irx   = collect(1:1:nx); irz = 1*ones(Int64, length(irx));
    shotv = MultiStepForward(irz, irx, srcs, fidMtx)
    return shotv
end

(vxx, vxz, vzx, vzz) = ReadPv(path_pv, 200);
p = (vzz+vxx); s = (vxz-vzx);
figure(); a = quantile(abs(p[:]),0.98); imshow(p, vmax=a, vmin=-a, cmap="PuOr");
figure(); a = quantile(abs(s[:]),0.98); imshow(s, vmax=a, vmin=-a, cmap="PuOr");

max_stars_repo_name: WneleiGao/ElasticWave


text:
abstract AbstractFElement{T <: AbstractMaterialStatus}

abstract AbstractElemStorage

@lazymod LinTrigMod "lin_trig.jl"
@lazymod LinQuadMod "lin_quad.jl"
@lazymod QuadTrigMod "quad_trig.jl"
@lazymod GradTrigMod "grad_trig.jl"
@lazymod GradTrigPrimalMod "grad_trig_prim.jl"

# Interface for a FE-element
createstorage() = error("Not Implemented")
createinterp() = error("Not Implemented")
Bmatrix() = error("Not Implemented")
creategps() = error("Not Implemented")
doftypes() = error("Not Implemented")
get_ndofs() = error("Not Implemented")
get_geotype() = error("Not Implemented")
get_ref_area() = error("Not Implemented")

getindex{T <: AbstractFElement}(elem::T, i0::Int) = getindex(elem.vertices, i0)
vertices(elem::AbstractFElement) = elem.vertices
gausspoints(elem::AbstractFElement) = elem.gps

function show{T <: AbstractFElement}(io::IO,elem::T)
    print(io, string(typeof(elem), ":", elem.vertices))
end

function stiffness(elem::AbstractFElement,
                  nodes::Vector{FENode2},
                  material::AbstractMaterial,
                  dof_vals)
    fill!(elem.storage.Ke, 0.0)
    for gp in elem.gps
        Be = Bmatrix(elem, gp, nodes)
        De = stiffness(material, gp)
        dV = weight(elem, gp, nodes)
        A_mul_B!(elem.storage.DeBe, De, Be) # DeBe = De * Be
        # Ke += B' * DeBe * dV
        BLAS.gemm!('T', 'N' ,dV, Be, elem.storage.DeBe, 1.0, elem.storage.Ke)
    end

    return elem.storage.Ke
end

# Extracts the values of the dofs for an element.
function get_field(elem::AbstractFElement, nodes::Vector{FENode2}, dof_vals::DofVals)
    u = elem.storage.u_field
    i = 1
    @inbounds for vert in elem.vertices
        for dof in nodes[vert].dofs
            u[i] = get_value(dof_vals, dof)
            i += 1
        end
    end
    return u
end

# Computes the internal forces for an element.
function intf(elem::AbstractFElement, mat::AbstractMaterial, nodes::Vector{FENode2}, dof_vals::DofVals)
    u = get_field(elem, nodes, dof_vals)
    ɛ = elem.storage.ɛ
    fill!(elem.storage.f_int, 0.0)
    for (i, gp) in enumerate(elem.gps)
        B = Bmatrix(elem, gp, nodes)
        A_mul_B!(ɛ, B, u)
        fill_from_start!(elem.temp_matstats[i].strain, ɛ)

        σ = stress(mat, ɛ, gp)
        fill_from_start!(elem.temp_matstats[i].stress, σ)

        dV = weight(elem, gp, nodes)

        # f_int += B' * σ * dV
        BLAS.gemv!('T', dV, B, σ, 1.0, elem.storage.f_int)
    end
    return elem.storage.f_int
end

function weight(elem::AbstractFElement, gp::GaussPoint2, nodes::Vector{FENode2})
    dN = dNmatrix(elem.interp, gp.local_coords)
    J = Jmatrix(elem.interp, elem.vertices, nodes, dN)
    return abs(det2x2(J)) * gp.weight
end

const ZERO_SCALAR = zeros(1)
const ZERO_VECTOR = zeros(3)
const ZERO_TENSOR = zeros(9)
get_field(elem::AbstractFElement, field::Type{AbstractScalar}, ::Int) = ZERO_SCALAR
get_field(elem::AbstractFElement, field::Type{AbstractVector}, ::Int) = ZERO_VECTOR
get_field(elem::AbstractFElement, field::Type{AbstractTensor}, ::Int) = ZERO_TENSOR

function get_cell_data{T <: AbstractField}(elem::AbstractFElement, field::Type{T})
    cellfield = zeros(get_ncomponents(field))
    for (i, gp) in enumerate(elem.gps)
        gpfield = get_field(elem, field, i)
        axpy!(getweight(gp)/get_ref_area(elem), gpfield, cellfield)
    end
    return cellfield
end

# Iterator for active dofs, too slow right now.
immutable ActiveDofsIterator{T <: AbstractFElement}
    nodes::Vector{FENode2}
    ele::T
end

get_field(elem::AbstractFElement, ::Type{Stress}, i::Int) = elem.matstats[i].stress
get_field(elem::AbstractFElement, ::Type{Strain}, i::Int) = elem.matstats[i].strain

function get_field(elem::AbstractFElement, ::Type{VonMises}, i::Int)
    σ = elem.matstats[i].stress
    m = (σ[1] + σ[2] + σ[3]) / 3
    return [sqrt(3/2) * sqrt((σ[1] - m)^2 + (σ[2] - m)^2 + (σ[3] - m)^2 +
                        2(σ[4]*σ[4]))]
end

function activedofs(element::AbstractFElement, nodes::Vector{FENode2})
    ActiveDofsIterator(nodes, element)
end

function Base.start(adi::ActiveDofsIterator)
    # We need to find first active dof because
    # we cant start the iterator loop unless
    # we actually have active dofs
    skip_inactive(adi, (1, 0, 0))
end

# Iterator for
function Base.done(adi::ActiveDofsIterator, state)
    state[1] > length(vertices(adi.ele))
end

function Base.next(adi::ActiveDofsIterator, state)
    v = state[1]
    i = state[2]
    n = state[3]

    node = adi.nodes[vertices(adi.ele)[v]]
    dof = get_dof(node, i)
    return((dof, n), skip_inactive(adi, state))
end

function skip_inactive(adi::ActiveDofsIterator, state)
    v = state[1]
    i = state[2] + 1
    n = state[3]

    node = adi.nodes[vertices(adi.ele)[v]]

    # Finish dofs for this node
    for i in i:length(get_dofs(node))
        n+= 1
        if isactive(get_dof(node, i))
            return (v, i, n)
        end
    end

    # Start looking for new active dof
    v+=1
    for v in v:length(vertices(adi.ele))
        node = adi.nodes[vertices(adi.ele)[v]]
        for i in 1:length(get_dofs(node))
            n += 1
            dof = get_dof(node, i)
            if isactive(dof)
                return (v, i, n)
            end
        end
    end
    return (v+1, i, n)
end


text:
module PermutedArrays

using StaticArrays
using Swizzles.WrapperArrays
using Swizzles.ShallowArrays
using Swizzles.ArrayifiedArrays
using Swizzles.StylishArrays
using Swizzles.ExtrudedArrays
using Swizzles
using Swizzles: SwizzledArray
using Swizzles: mask, jointuple, combinetuple, masktuple, ziptuple, nziptuple

using Base: Slice
using Base.Broadcast: BroadcastStyle, Broadcasted, Style, AbstractArrayStyle, DefaultArrayStyle
using Base.Broadcast: broadcast_shape, check_broadcast_shape, result_style, axistype, broadcasted
using Base.Iterators: repeated, flatten, take

export PermuteStyle, Permute, PermutationMismatch

"""
    isperm(v::AbstractVector) -> Bool

Return `true` if `v` is a valid permutation of `axes(v, 1)`.
"""
Base.isperm(v::AbstractVector) = Base.isperm(v, axes(v, 1))

"""
    isperm(v, u) -> Bool

Return `true` if `v` is a valid permutation of `u`.
"""
function Base.isperm(v, u)
    return sort(v) == sort(u)
end
function Base.isperm(v::UnitRange{T}, u::UnitRange{T}) where {T<:Integer}
    return v == u
end
function Base.isperm(v::UnitRange{T}, u::AbstractVector{T}) where {T<:Integer}
    return isperm(u, v)
end
function Base.isperm(v::AbstractVector{T}, u::UnitRange{T}) where {T<:Integer}
    n = length(u)
    used = falses(n)
    for x in v
        i = x - u.start
        (0 < x <= n) && (used[x] ⊻= true) || return false
    end
    true
end
function Base.isperm(v::Base.OneTo{T}, u::Base.OneTo{T}) where {T<:Integer}
    return v == u
end
function Base.isperm(v::Base.OneTo{T}, u::AbstractVector{T}) where {T<:Integer}
    return isperm(u, v)
end
function Base.isperm(v::AbstractVector{T}, u::Base.OneTo{T}) where {T<:Integer}
    n = length(u)
    used = falses(n)
    for x in v
        i = x - 1
        (0 < x <= n) && (used[x] ⊻= true) || return false
    end
    true
end

Base.invperm(v::UnitRange{<:Integer}) = v
Base.invperm(v::Base.OneTo{<:Integer}) = v

"""
    PermutationMismatch([msg])

The objects called do not have matching permutations. Optional argument `msg` is
a descriptive error string.
"""
struct PermutationMismatch <: Exception
    msg::AbstractString
end
PermutationMismatch() = PermutationMismatch("")

"""
permute an array so that
`arr[i...] = parent(arr)[map(getindex, perms(arr), i)]
interesting observation: perms[i] must be a permutation of axes(perms[i]).
"""
struct PermutedArray{T, N, Perms <: Tuple{Vararg{AbstractVector, N}}, IPerms <: Tuple{Vararg{AbstractVector, N}}, Arg <: AbstractArray} <: StylishArray{T, N}
    perms::Perms
    iperms::IPerms
    arg::Arg
end

PermutedArray(perms, iperms, arg) = PermutedArray{eltype(arg), ndims(arg), typeof(perms), typeof(iperms), typeof(arg)}(perms, iperms, arg)

Base.parent(arr::PermutedArray) = arr.arg
WrapperArrays.iswrapper(arr::PermutedArray) = true
WrapperArrays.adopt(arr::PermutedArray, arg) = PermutedArray(arr.perms, arr.iperms, arg)

Base.Broadcast.BroadcastStyle(::Type{<:PermutedArray{<:Any, <:Any, <:Any, <:Any, Arg}}) where {Arg} = PermuteStyle(BroadcastStyle(Arg))

Base.IndexStyle(arr::PermutedArray) = IndexCartesian()
Base.IndexStyle(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.OneTo, N}}}) where {T, N} = IndexStyle(arr.arg)
Base.IndexStyle(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.Slice{<:Base.OneTo}, N}}}) where {T, N} = IndexStyle(arr.arg)
Base.size(arr::PermutedArray) = size(arr.arg)
Base.size(arr::PermutedArray, d::Int) = size(arr.arg, d)
Base.axes(arr::PermutedArray) = axes(arr.arg)
Base.axes(arr::PermutedArray, d::Int) = axes(arr.arg, d)
(Base.getindex(arr::PermutedArray{T, 0, Tuple{}})::T) where {T} = arr.arg[]
function Base.getindex(arr::PermutedArray{T, N}, i::Vararg{Any, N})::T where {T, N}
    arr.arg[ntuple(n->arr.perms[n][i[n]], Val(N))...]
end
(Base.getindex(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.OneTo, N}}}, i...)::T) where {T, N} = arr.arg[i...]
(Base.getindex(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.Slice{<:Base.OneTo}, N}}}, i...)::T) where {T, N} = arr.arg[i...]
(Base.getindex(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.OneTo, N}}}, i::Vararg{Any, N})::T) where {T, N} = arr.arg[i...]
(Base.getindex(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.Slice{<:Base.OneTo}, N}}}, i::Vararg{Any, N})::T) where {T, N} = arr.arg[i...]
(Base.getindex(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.OneTo, N}}}, i::Integer)::T) where {T, N} = arr.arg[i]
(Base.getindex(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.Slice{<:Base.OneTo}, N}}}, i::Integer)::T) where {T, N} = arr.arg[i]
(Base.setindex!(arr::PermutedArray{T, 0, Tuple{}}, x)::T) where {T} = arr.arg[] = x
function Base.setindex!(arr::PermutedArray{T, N}, x, i::Vararg{Any, N})::T where {T, N}
    arr.arg[ntuple(n->arr.iperms[n][i[n]], Val(N))...] = x
end
(Base.setindex!(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.OneTo, N}}}, x, i::Vararg{Any, N})::T) where {T, N} = arr.arg[i...] = x
(Base.setindex!(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.Slice{<:Base.OneTo}, N}}}, x, i::Vararg{Any, N})::T) where {T, N} = arr.arg[i...] = x
(Base.setindex!(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.OneTo, N}}}, x, i...)::T) where {T, N} = arr.arg[i...] = x
(Base.setindex!(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.Slice{<:Base.OneTo}, N}}}, x, i...)::T) where {T, N} = arr.arg[i...] = x
(Base.setindex!(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.OneTo, N}}}, x, i::Integer)::T) where {T, N} = arr.arg[i] = x
(Base.setindex!(arr::PermutedArray{T, N, <:Tuple{Vararg{<:Base.Slice{<:Base.OneTo}, N}}}, x, i::Integer)::T) where {T, N} = arr.arg[i] = x

struct PermuteStyle{S<:BroadcastStyle} <: BroadcastStyle end
PermuteStyle(style::S) where {S <: BroadcastStyle} = PermuteStyle{S}()
PermuteStyle(style::S) where {S <: PermuteStyle} = style
Base.Broadcast.BroadcastStyle(::PermuteStyle{T}, ::S) where {T, S<:AbstractArrayStyle} = PermuteStyle(result_style(T(), S()))
Base.Broadcast.BroadcastStyle(::PermuteStyle{T}, ::S) where {T, S<:DefaultArrayStyle} = PermuteStyle(result_style(T(), S()))
Base.Broadcast.BroadcastStyle(::PermuteStyle{T}, ::S) where {T, S<:Style{Tuple}} = PermuteStyle(result_style(T(), S()))
Base.Broadcast.BroadcastStyle(::PermuteStyle{T}, ::PermuteStyle{S}) where {T, S<:BroadcastStyle} = PermuteStyle(result_style(T(), S()))

Base.@propagate_inbounds function Base.similar(bc::Broadcasted{PermuteStyle{S}}, args...) where {S}
    src = repermute(bc)
    return PermutedArray(src.perms, src.iperms, similar(convert(Broadcasted{S}, bc), args...))
end

Base.@propagate_inbounds function Base.copy(bc::Broadcasted{PermuteStyle{S}}) where {S <: Base.Broadcast.AbstractArrayStyle}
    res = copyto!(similar(bc), bc)
    return ndims(res) == 0 ? res[] : res
end

Base.@propagate_inbounds function Base.copyto!(dst::AbstractArray, src::Broadcasted{PermuteStyle{S}}) where {S}
    pdst = repermute(dst)
    psrc = repermute(src)
    pdst.perms == psrc.perms || PermutationMismatch("TODO")
    copyto!(pdst.arg, psrc.arg)
    return dst
end

@inline function repermute(src)
    perms = map(Slice, axes(src))
    iperms = map(invperm, perms)
    return PermutedArray(perms, iperms, arrayify(src))
end

@inline function repermute(src::PermutedArray)
    arg = repermute(src.arg)
    return PermutedArray(map(getindex, src.perms, arg.perms), map(getindex, src.iperms, arg.iperms), arg.arg)
end

@inline function repermute(src::Broadcasted{PermuteStyle{S}}) where {S}
    return repermute(convert(Broadcasted{S}, src))
end
@inline function repermute(src::Broadcasted{S}) where {S}
    args = map(repermute, src.args)
    combineables = map(arg->ziptuple(arg.perms, arg.iperms, keeps(arg)), args)
    (perms, iperms, _) = nziptuple(Val(3), combinetuple(result_perm, combineables...)...)
    args = map(parent, args)
    return PermutedArray(perms, iperms, ArrayifiedArray(Broadcasted{S}(src.f, args, src.axes)))
end

function result_perm((p1, ip1, k1), (p2, ip2, k2))
    if !kept(k2)
        return (permtype(p1, p2), permtype(ip1, ip2), k1 | k2)
    elseif !kept(k1)
        return (permtype(p2, p1), permtype(ip2, ip1), k1 | k2)
    else
        p1 == p2 || throw(PermutationMismatch("Conflicting permutations"))
        return (permtype(p1, p2), permtype(ip1, ip2), k1 | k2)
    end
end
@inline result_perm((p1, ip1, k1)::Tuple{Any, Keep}, (p2, ip2, k2)::Tuple{Any, Extrude}) = (p1, ip1, Keep())
@inline result_perm((p1, ip1, k1)::Tuple{Any, Extrude}, (p2, ip2, k2)::Tuple{Any, Keep}) = (p2, ip2, Keep())
@inline function result_perm((p1, ip1, k1)::Tuple{Any, Keep}, (p2, ip2, k2)::Tuple{Any, Keep})
    p1 == p2 || throw(PermutationMismatch("Conflicting permutations"))
    (p1, ip1, Keep())
end

@inline permtype(p1::P, p2::P) where {P} = p1
@inline permtype(p1::AbstractVector{Int}, p2::AbstractVector{Int}) = Vector{Int}(p1)
@inline permtype(p1::Base.OneTo, p2::Base.OneTo) = Base.OneTo(convert(promote_type(eltype(p1), eltype(p2)), p1[end]))

@inline Swizzles.childstyle(Arr::Type{<:SwizzledArray}, ::PermuteStyle{S}) where {S} = PermuteStyle(Swizzles.childstyle(Arr, S()))

function repermute(arr::SwizzledArray{T}) where {T}
    arg = repermute(arr.arg)::PermutedArray
    init = repermute(arr.init)::PermutedArray
    arr_perms = masktuple(d->Base.OneTo(1), d->arg.perms[d], Val(mask(arr)))
    arr_iperms = masktuple(d->Base.OneTo(1), d->arg.iperms[d], Val(mask(arr)))
    arg_keeps = keeps(arg)
    arr_keeps = masktuple(d->Extrude(), d->arg_keeps[d], Val(mask(arr)))
    init_perms = init.perms
    init_iperms = init.iperms
    init_keeps = keeps(init)
    combineables = (ziptuple(arr_perms, arr_iperms, arr_keeps), ziptuple(init_perms, init_iperms, init_keeps))
    (perms, iperms, _) = nziptuple(Val(3), combinetuple(result_perm, combineables...)...)
    return PermutedArray(perms, iperms, Swizzle(arr.op, Val(mask(arr)))(parent(init), parent(arg)))
end

repermute(arr::ArrayifiedArray) = repermute(arr.arg)

struct Permute{_Perms, _IPerms} <: Swizzles.Intercept
    _perms::_Perms
    _iperms::_IPerms
    @inline function Permute(_perms...)
        _iperms = map(invperm, _perms)
        return new{typeof(_perms), typeof(_iperms)}(_perms, _iperms)
    end
end

@inline function(ctr::Permute{<:Tuple{Vararg{Any, _N}}, <:Tuple{Vararg{Any, _N}}})(arg::Arg) where {_N, Arg <: AbstractArray}
    if @generated
        perms = ntuple(n -> n <= _N ? :(ctr._perms[$n]) : :(axes(arg, $n)), ndims(arg))
        iperms = ntuple(n -> n <= _N ? :(ctr._iperms[$n]) : :(axes(arg, $n)), ndims(arg))
        return :(return PermutedArray(($(perms...),), ($(iperms...),), arg))
    else
        perms = ntuple(n -> n <= _N ? ctr._perms[n] : axes(arg, n), ndims(arg))
        iperms = ntuple(n -> n <= _N ? ctr._iperms[n] : axes(arg, n), ndims(arg))
        return PermutedArray(perms, iperms, arg)
    end
end

end


text:
## Assignments.jl --- for secondary structure assignments/classifications

"""
    abstract type AssignmentScheme

Supertype of all secondary structure assignment schemes. Plans are to implement
DSSP, STRIDE, KAKSI, (three-state) HEO schemes.
"""
abstract type AssignmentScheme end

struct DSSP <: AssignmentScheme end
struct STRIDE <: AssignmentScheme end
struct KAKSI <: AssignmentScheme end
struct HEO <: AssignmentScheme end

### TODO: declare enums for secondary structures assigned by each scheme
### TODO: abstraction for creating and applying assignment criteria
### TODO: calculation of criteria specifics (in respective files)


text:
using JuMP, Ipopt, NLsolve, BenchmarkTools

# Ipopt approach
numSectors = 25
u(c) = log(c)
modTrial = Model(Ipopt.Optimizer)
set_silent(modTrial)
@variable(modTrial, c[1:numSectors] >= 0.01)
@NLobjective(modTrial, Min, sum((u(c[i]))^2 for i in 1:numSectors))
@btime JuMP.optimize!(modTrial)
println(value.(c))

# NLsolve approach
function f!(F, x)
    for i in 1:numSectors
        F[i] = u(x[i])
    end
end
initial_x = 1.5*ones(numSectors)
@btime nlsolve(f!, initial_x)

#NLsolve seems faster

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
using LogicCircuits: is⋁gate, is⋀gate

function test_pc_equals(c1, c2)
    @test num_nodes(c1) == num_nodes(c2)
    @test num_edges(c1) == num_edges(c2)
    for (n1, n2) in zip(linearize(c1), linearize(c2))
        if is⋁gate(n1)
            @test is⋁gate(n2)
            @test all(params(n1) ≈ params(n2))
        elseif is⋀gate(n1)
            @test is⋀gate(n2)
        else
            @test isliteralgate(n1)
            @test isliteralgate(n2)
            @test literal(n1) == literal(n2)
        end
    end
end


text:
#
# Calculates the SVD of a data matrix containing missing entries.
#
# Uses the iterative algorithm of Hastie et al. 1999
#

# Calculate the rank-k SVD approximation to a matrix given the
# full SVD.
function approximate(u, d, v, k::Int64)
  u[:, 1:k] * diagm(d[1:k]) * v[1:k, :]
end

# Test code:
#
# srand(1)
# M = rand(3, 3)
# (u, d, v) = svd(M)
# norm(approximate(u, d, v, 1) - M)
# norm(approximate(u, d, v, 2) - M)
# norm(approximate(u, d, v, 3) - M)

# Impute a missing entries using current approximation.
function impute(m::Matrix{Float64}, missing_entries, u, d, v, k::Int64)
  approximate_m = approximate(u, d, v, k)

  for index_pair in missing_entries
    m[index_pair[1], index_pair[2]] = approximate_m[index_pair[1], index_pair[2]]
  end

  m
end

# Should be done with a proper N-dimensional Int64 array.
function ind_na{T}(df::DataMatrix{T})
  indices = {}
  for i = 1:nrow(df)
    for j = 1:ncol(df)
      if isna(df[i, j])
        push!(indices, [i, j])
      end
    end
  end
  indices
end

# Kind of a nutty method without DataMatrix.
function mean{T}(df::DataMatrix{T})
  mu = 0.0
  n = 0
  for i = 1:nrow(df)
    for j = 1:ncol(df)
      if !isna(df[i, j])
        mu += df[i, j]
        n += 1
      end
    end
  end
  mu / n
end

# This will crash if a row is missing all entries.
function rowmeans{T}(df::DataMatrix{T})
  mus = zeros(nrow(df))
  for i = 1:nrow(df)
    mu = 0.0
    n = 0
    for j = 1:ncol(df)
      if !isna(df[i, j])
        mu += df[i, j]
        n += 1
      end
    end
    mus[i] = mu / n
  end
  mus
end

# Must select rank k of SVD to use.
# TODO: Default to failure in the face of NA's
function svd{T}(D::DataMatrix{T}, k::Int)
  df = deepcopy(D)

  print_trace = false

  tolerance = 10e-4

  # Cache the dimensions of the matrix.
  n = size(df, 1)
  p = size(df, 2)

  # Estimate missingness and print a message.
  missing_entries = ind_na(df)
  missingness = length(missing_entries) / (nrow(df) * ncol(df))
  if print_trace
    println("Matrix is missing $(missingness * 100)% of entries")
  end

  # Initial imputation uses row means.
  global_mu = mean(df)
  mu_i = rowmeans(df)

  for i = 1:n
    for j = 1:p
      if isna(df[i, j])
        if isna(mu_i[i])
          df[i, j] = global_mu
        else
          df[i, j] = mu_i[i]
        end
      end
    end
  end

  # Make a matrix out of the dataframe.
  tmp = zeros(n, p)
  for i = 1:n
    for j = 1:p
      tmp[i, j] = df[i, j]
    end
  end
  df = tmp

  # Count iterations of proper imputation method.
  i = 0

  # Keep track of approximate matrices.
  previous_df = copy(df)
  current_df = copy(df)

  # Keep track of Frobenius norm of changes in imputed matrix.
  change = Inf

  # Iterate until imputation stops changing up to a tolerance of 10e-6.
  while change > tolerance
    if print_trace
      println("Iteration $i")
      println("Change $change")
    end

    # Impute missing entries using current SVD.
    previous_df = copy(current_df)
    u, d, v = svd(current_df)
    current_df = impute(current_df, missing_entries, u, d, v, k)

    # Compute the change in the matrix across iterations.
    change = norm(previous_df - current_df) / norm(df)

    # Increment the iteration counter.
    i = i + 1
  end

  # Tell the user how many iterations were required to impute matrix.
  if print_trace
    println("Tolerance achieved after $i iterations")
  end

  # Return both df and the SVD of df with all entries imputed.
  u, d, v = svd(current_df)

  # Only return the SVD entries, not the imputation
  return (u[:, 1:k], d[1:k], v[1:k, :])
end

svd{T}(D::DataMatrix{T}) = svd(D, min(nrow(D), ncol(D)))

function eig{T}(D::DataMatrix{T})
  u, d, v = svd(D)
  eig(u * diagm(d) * v)
end


text:
"Holds the tableau of an Specialised Partitioned Additive Runge-Kutta method for Variational systems."
struct TableauVSPARKsecondary{DT <: Number} <: AbstractTableau{DT}
    name::Symbol
    o::Int
    s::Int
    r::Int
    ρ::Int

    q::CoefficientsSPARK{DT}
    p::CoefficientsSPARK{DT}

    q̃::CoefficientsSPARK{DT}
    p̃::CoefficientsSPARK{DT}

    ω::Matrix{DT}
    d::Vector{DT}

    function TableauVSPARKsecondary(name::Symbol, o::Int, s::Int, r::Int,
                        q::CoefficientsSPARK{DT}, p::CoefficientsSPARK{DT},
                        q̃::CoefficientsSPARK{DT}, p̃::CoefficientsSPARK{DT},
                        ω::Matrix{DT}, d=DT[]) where {DT}

        @assert s > 0 "Number of stages s must be > 0"
        @assert r > 0 "Number of stages r must be > 0"

        @assert s==q.s==p.s==q̃.σ==p̃.σ
        @assert r==q.σ==p.σ==q̃.s==p̃.s
        @assert size(ω,1)==r
        @assert size(ω,2)==r+1

        @assert length(d)==0 || length(d)==r

        new{DT}(name, o, s, r, 0, q, p, q̃, p̃, ω, d)
    end
end

"Parameters for right-hand side function of Specialised Partitioned Additive Runge-Kutta methods for Variational systems."
const ParametersVSPARKsecondary = AbstractParametersSPARK{:vspark_secondary}

@doc raw"""
Specialised Partitioned Additive Runge-Kutta integrator for degenerate
variational systems with projection on secondary constraint.

This integrator solves the following system of equations for the internal stages,
```math
\begin{aligned}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a^1_{ij} V_{n,j} + h \sum \limits_{j=1}^{\sigma} a^2_{ij} \tilde{\Lambda}_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \bar{a}^1_{ij} F^1_{n,j} + h \sum \limits_{j=1}^{s} \bar{a}^2_{ij} F^2_{n,j} + h \sum \limits_{j=1}^{\sigma} \bar{a}^3_{ij} \tilde{F}^3_{n,j} , & i &= 1, ..., s , \\
0 &= \Phi_{n,i} , & i &= 1, ..., s ,
\end{aligned}
```
the projective stages
```math
\begin{aligned}
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \alpha^1_{ij} \tilde{V}_{n,j} + h \sum \limits_{j=1}^{\sigma} \alpha^2_{ij} \tilde{\Lambda}_{n,j} , & i &= 1, ..., \sigma , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \bar{\alpha}^1_{ij} F^1_{n,j} + h \sum \limits_{j=1}^{s} \bar{\alpha}^2_{ij} \tilde{F}^2_{n,j} + h \sum \limits_{j=1}^{\sigma} \bar{\alpha}^3_{ij} \tilde{F}^3_{n,j} , & i &= 1, ..., \sigma , \\
0 &= \tilde{\Phi}_{n,i} , & i &= 1, ..., \sigma , \\
0 &= \sum \limits_{j=1}^{\sigma} \omega_{ij} \tilde{\Psi}_{n,i} , & i &= 1, ..., \sigma-1 ,
\end{aligned}
```
and update rule
```math
\begin{aligned}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b^1_{i} V_{n,i} + h \sum \limits_{i=1}^{\sigma} b^2_{i} \tilde{\Lambda}_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b^1_{i} F^1_{n,i} + h \sum \limits_{i=1}^{s} b^2_{i} F^2_{n,i} + h \sum \limits_{i=1}^{\sigma} b^3_{i} \tilde{F}^3_{n,i} , \\
0 &= \phi (q_{n+1}, p_{n+1}) ,
\end{aligned}
```
with definitions
```math
\begin{aligned}
F^1_{n,i} &= \nabla H (Q_{n,i}) , & i &= 1, ..., s , \\
F^2_{n,i} &= \nabla \vartheta (Q_{n,i}) \cdot V_{n,i} , & i &= 1, ..., s , \\
\tilde{F}^3_{n,i} &= \nabla \phi (\tilde{Q}_{n,i}, \tilde{P}_{n,i}) \cdot \tilde{\Lambda}_{n,i} , & i &= 1, ..., \sigma , \\
\Phi_{n,i} &= \phi (Q_{n,i}, P_{n,i}) = P_{n,i} - \vartheta (Q_{n,i}) , & i &= 1, ..., s , \\
\tilde{\Phi}_{n,i} &= \phi (\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = \tilde{P}_{n,i} - \vartheta (\tilde{Q}_{n,i}) , & i &= 1, ..., \sigma , \\
\tilde{\Psi}_{n,i} &= \psi (\tilde{Q}_{n,i}, \tilde{V}_{n,i}, \tilde{P}_{n,i}, \tilde{F}_{n,i}) = \tilde{F}_{n,i} - \tilde{V}_{n,i} \cdot \nabla \vartheta (\tilde{Q}_{n,i}) , & i &= 1, ..., \sigma ,
\end{aligned}
```
so that
```math
\begin{aligned}
F^1_{n,i} + F^2_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
\phi (Q_{n,i}, P_{n,i}) &= P_{n,i} - \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
\psi (\tilde{Q}_{n,i}, \tilde{V}_{n,i}, \tilde{P}_{n,i}, \tilde{F}_{n,i}) &= \tilde{F}_{n,i} - \tilde{V}_{n,i} \cdot \nabla \frac{\partial L}{\partial v} (\tilde{Q}_{n,i}, \tilde{V}_{n,i}) , & i &= 1, ..., \sigma .
\end{aligned}
```
"""
struct IntegratorVSPARKsecondary{DT, TT, D, S, R,
                                         PT <: ParametersVSPARKsecondary{DT,TT,D,S,R},
                                         ST <: NonlinearSolver{DT},
                                         IT <: InitialGuessIODE{TT}} <: AbstractIntegratorVSPARK{DT,TT,D,S,R}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorVSPARKsecondary(params::ParametersVSPARKsecondary{DT,TT,D,S,R}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,R,ST,IT}
        new{DT, TT, D, S, R, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorVSPARKsecondary{DT,D}(equations::NamedTuple, tableau::TableauVSPARKsecondary{TT}, Δt::TT) where {DT,TT,D,ST}
        # get number of stages
        S = tableau.s
        R = tableau.r

        N = 4*D*R

        if isdefined(tableau, :d) && length(tableau.d) > 0
            N += D
        end

        # create params
        params = ParametersVSPARKsecondary{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, N, params, caches)

        # create initial guess
        iguess = InitialGuessIODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt)

        # create integrator
        IntegratorVSPARKsecondary(params, solver, iguess, caches)
    end

    function IntegratorVSPARKsecondary(equation::LDAE{DT}, tableau::TableauVSPARKsecondary, Δt; kwargs...) where {DT}
        IntegratorVSPARKsecondary{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

GeometricBase.nconstraints(::IntegratorVSPARKsecondary{DT,TT,D}) where {DT,TT,D} = D

function initial_guess!(int::IntegratorVSPARKsecondary{DT}, sol::AtomicSolutionPDAE{DT},
                        cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT}
    for i in 1:pstages(int)
        evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                              sol.q, sol.p, sol.v, sol.f,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in 1:ndims(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - sol.p[k])/timestep(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+3] = cache.ṽ[k]
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+4] = 0
        end
    end

    if isdefined(tableau(int), :d) && length(tableau(int).d) > 0
        for k in 1:ndims(int)
            int.solver.x[4*ndims(int)*pstages(int)+k] = 0
        end
    end
end

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S,R},
                                        params::ParametersVSPARKsecondary{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local t::TT

    for i in 1:R
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[4*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[4*(D*(i-1)+k-1)+2]
            cache.Vp[i][k] = x[4*(D*(i-1)+k-1)+3]
            cache.Λp[i][k] = x[4*(D*(i-1)+k-1)+4]

            # compute Q and P
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        t = params.t + params.Δt * params.tab.p̃.c[i]
        params.equs[:f](t, cache.Qp[i], cache.Vp[i], cache.Fp[i])
        params.equs[:g](t, cache.Qp[i], cache.Vp[i], cache.Gp[i])
        params.equs[:g](t, cache.Qp[i], cache.Λp[i], cache.G̅p[i])

        cache.Hp[i] .= cache.Fp[i] .+ cache.Gp[i]

        params.equs[:ϕ](t, cache.Qp[i], cache.Pp[i], cache.Φp[i])
        params.equs[:ψ](t, cache.Qp[i], cache.Pp[i], cache.Vp[i], cache.Hp[i], cache.Ψp[i])
    end

    if isdefined(params.tab, :d) && length(params.tab.d) > 0
        for k in 1:D
            cache.μ[k] = x[4*D*R+k]
        end
    end

    for i in 1:S
        # compute Q
        for k in 1:D
            cache.Yi[i][k] = 0
            for j in 1:R
                cache.Yi[i][k] += params.tab.q.a[1][i,j] * cache.Vp[j][k]
                cache.Yi[i][k] += params.tab.q.a[2][i,j] * cache.Λp[j][k]
            end
            cache.Qi[i][k] = params.q[k] + params.Δt * cache.Yi[i][k]
        end

        # compute f(X)
        t = params.t + params.Δt * params.tab.p.c[i]
        params.equs[:f](t, cache.Qi[i], cache.Vi[i], cache.Fi[i])
    end

    # compute q and p
    cache.q̃ .= params.q
    cache.p̃ .= params.p
    for i in 1:S
        cache.p̃ .+= params.Δt .* params.tab.p.b[1][i] .* cache.Fi[i]
    end
    for i in 1:R
        cache.q̃ .+= params.Δt .* params.tab.q.b[1][i] .* cache.Vp[i]
        cache.q̃ .+= params.Δt .* params.tab.q.b[2][i] .* cache.Λp[i]
        cache.p̃ .+= params.Δt .* params.tab.p.b[2][i] .* cache.Gp[i]
        cache.p̃ .+= params.Δt .* params.tab.p.b[3][i] .* cache.G̅p[i]
    end

    # compute ϕ(q,p)
    t = params.t + params.Δt
    params.equs[:ϕ](t, cache.q̃, cache.p̃, cache.ϕ̃)
end

"Compute stages of specialised partitioned additive Runge-Kutta methods for variational systems."
function Integrators.function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersVSPARKsecondary{DT,TT,D,S,R},
                                      caches::CacheDict) where {ST,DT,TT,D,S,R}

    # get cache for internal stages
    cache = caches[ST]

    compute_stages!(y, cache, params)

    # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ, ωΨ]
    for i in 1:R
        for k in 1:D
            b[4*(D*(i-1)+k-1)+1] = - cache.Yp[i][k]
            b[4*(D*(i-1)+k-1)+2] = - cache.Zp[i][k]
            b[4*(D*(i-1)+k-1)+3] = - cache.Φp[i][k]
            b[4*(D*(i-1)+k-1)+4] = 0
            for j in 1:S
                b[4*(D*(i-1)+k-1)+2] += params.tab.p̃.a[1][i,j] * cache.Fi[j][k]
            end
            for j in 1:R
                b[4*(D*(i-1)+k-1)+1] += params.tab.q̃.a[1][i,j] * cache.Vp[j][k]
                b[4*(D*(i-1)+k-1)+1] += params.tab.q̃.a[2][i,j] * cache.Λp[j][k]
                b[4*(D*(i-1)+k-1)+2] += params.tab.p̃.a[2][i,j] * cache.Gp[j][k]
                b[4*(D*(i-1)+k-1)+2] += params.tab.p̃.a[3][i,j] * cache.G̅p[j][k]
            end
            for j in 1:R
                b[4*(D*(i-1)+k-1)+4] -= params.tab.ω[i,j] * cache.Ψp[j][k]
            end
            b[4*(D*(i-1)+k-1)+4] -= params.tab.ω[i,R+1] * cache.ϕ̃[k]
        end
    end

    if isdefined(params.tab, :d) && length(params.tab.d) > 0
        for i in 1:R
            for k in 1:D
                b[4*(D*(i-1)+k-1)+2] -= cache.μ[k] * params.tab.d[i] / params.tab.p.b[2][i]
            end
        end

        for k in 1:D
            b[4*D*R+k] = 0
            for i in 1:R
                b[4*D*R+k] -= cache.Vp[i][k] * params.tab.d[i]
            end
        end
    end
end

function update_solution!(int::IntegratorVSPARKsecondary{DT,TT}, sol::AtomicSolutionPDAE{DT,TT},
                          cache::IntegratorCacheSPARK{DT}=int.caches[DT]) where {DT,TT}
    # compute final update
    update_solution!(sol.p, sol.p̃, cache.Fi, int.params.tab.p.b[1], timestep(int))

    # compute projection
    update_solution!(sol.q, sol.q̃, cache.Vp, int.params.tab.q.b[1], timestep(int))
    update_solution!(sol.q, sol.q̃, cache.Λp, int.params.tab.q.b[2], timestep(int))
    update_solution!(sol.p, sol.p̃, cache.Gp, int.params.tab.p.b[2], timestep(int))
    update_solution!(sol.p, sol.p̃, cache.G̅p, int.params.tab.p.b[3], timestep(int))
end


text:
"""
    FileStack{X,K}

    FileStack{X,K}(filename, types, sizes, eachchunk, haschunks, write)

A wrapper object that holds file pointer and size/chunking
metadata for a multi-layered stack stored in a single file,
typically netcdf or hdf5.

`X` is a backend singleton like `NCDfile`, and `K` is a tuple
of `Symbol` keys.
"""
struct FileStack{X,K,F<:AbstractString,T,S,EC,HC}
    filename::F
    types::T
    sizes::S
    eachchunk::EC
    haschunks::HC
    write::Bool
end
function FileStack{X,K}(
    filename::F, types::T, sizes::S, eachchunk::EC, haschunks::HC, write::Bool
   ) where {X,K,F,T,S,EC,HC}
    FileStack{X,K,F,T,S,EC,HC}(filename, types, sizes, eachchunk, haschunks, write)
end

# FileStack has `X` parameter that is not recoverable from fields.
ConstructionBase.constructorof(::Type{<:FileStack{X,K}}) where {X,K} = FileStack{X,K}

filename(fs::FileStack) = fs.filename
Base.keys(fs::FileStack{<:Any,K}) where K = K
Base.values(fs::FileStack{<:Any}) = (fs[k] for k in keys(fs))
# Indexing FileStack returns a FileArray,
# referencing a specific key in the same file.
function Base.getindex(fs::FileStack{X,K}, key::Symbol) where {X,K}
    is = NamedTuple{K}(ntuple(identity, length(K)))
    i = is[key]
    size = fs.sizes[i]
    eachchunk = fs.eachchunk[i]
    haschunks = fs.haschunks[i]
    T = fs.types[i]
    N = length(size)
    A = FileArray{X,T,N}(filename(fs), size, key, eachchunk, haschunks, fs.write)
end


text:
using Documenter
using ConstrainedDynamics

makedocs(;
    modules = [ConstrainedDynamics],
    format = Documenter.HTML(
        canonical = "https://janbruedigam.github.io/ConstrainedDynamics.jl/stable/",
        assets = ["assets/favicon.ico"],
    ),
    pages = [
        "Home" => "index.md",
        "Getting started" => [
            "From URDF" => "gettingstarted/urdf.md",
            "From Code" => "gettingstarted/code.md"
        ],
        "Library" => [
            "Bodies and Origin" => "library/body.md",
            "Constraints" => "library/constraint.md",
            "Mechanism" => "library/mechanism.md",
            "Interface" => "library/interface.md",
            "Simulation" => "library/simulation.md",
        ]
    ],
    sitename = "ConstrainedDynamics.jl"
)

deploydocs(; repo = "github.com/janbruedigam/ConstrainedDynamics.jl")


text:
struct NoiseProblem{N<:AbstractNoiseProcess,T} <: AbstractNoiseProblem
  noise::N
  tspan::T
  seed::UInt64
end

@add_kwonly function NoiseProblem(noise,tspan;seed=UInt64(0))
  _tspan = promote_tspan(tspan)
  NoiseProblem{typeof(noise),typeof(_tspan)}(noise,_tspan,seed)
end


text:
struct Constant <: Expression
    ref::MyRef{Float64}
    Constant(x::R) where R <: Real  = new(MyRef(x))
end
struct Variable <: Expression
    index::Int
    ref::MyRef{Float64}
    Variable(n::Int) = new(n,MyRef(0.))
end
struct Parameter <: Expression
    index::Int
    ref::MyRef{Float64}
    Parameter(n::Int) = new(n,MyRef(0.))
end

struct Expression1{F <: Function ,E <: Expression} <: Expression
    e1::E
    ref::MyRef{Float64}
    Expression1(f::F,e1::E) where {F,E} = new{F,E}(e1,MyRef(0.))
end
struct Expression2{F <: Function,E1,E2} <: Expression
    e1::E1
    e2::E2
    ref::MyRef{Float64}
    Expression2(f::F,e1::E1,e2::E2) where {F,E1,E2} = new{F,E1,E2}(e1,e2,MyRef(0.))
end
struct ExpressionSum{E <: Expression,I} <: Expression
    inner::I
    es::Vector{E}
    ref::MyRef{Float64}
    ExpressionSum(es) = new{eltype(es),Nothing}(nothing,es,MyRef(0.))
    ExpressionSum(inner,es) = new{eltype(es),typeof(inner)}(inner,es,ref(inner))
end
struct Equality{E<:Expression} e::E end
struct Inequality{E<:Expression} e::E end

@inline function (e::ExpressionSum{E,I})(x,p=nothing) where {E,I}
    inner(e)(x,p)
    @simd for i in eachindex(e.es)
        @inbounds addrefval(e,e.es[i](x,p))
    end
    return refval(e)
end
@inline function (e::ExpressionSum{E,Nothing})(x,p=nothing) where E
    setrefval(e,.0)
    @simd for i in eachindex(e.es)
        @inbounds addrefval(e,e.es[i](x,p))
    end
    return refval(e)
end
@inline (e::Variable)(x,p=nothing)  = setrefval(e,getindex(x,index(e)))
@inline (e::Parameter)(x,p=nothing)  = setrefval(e,getindex(p,index(e)))
@inline (e::Constant)(x,p=nothing)  = refval(e)

for (T1,T2) in [(Expression,Expression),(Expression,Real),(Real,Expression)]
    @eval begin
        ==(e1::T1,e2::T2) where {T1 <: $T1, T2 <: $T2} = Equality(e1-e2)
        >=(e1::T1,e2::T2) where {T1 <: $T1, T2 <: $T2} = Inequality(e1-e2)
        <=(e1::T1,e2::T2) where {T1 <: $T1, T2 <: $T2} = Inequality(e2-e1)
    end
end


text:
@doc raw"""
    gausslobatto(n::Integer) -> Tuple{Vector{Float64},Vector{Float64}}

Return nodes and weights of [Gauss-Lobatto quadrature](https://mathworld.wolfram.com/LobattoQuadrature.html).

```math
\int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i)
```

# Examples
```jldoctest
julia> x, w = gausslobatto(4);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 2/5
true
```

Note that the both ends of nodes are fixed at -1 and 1.

```jldoctest
julia> x, w = gausslobatto(4);

julia> x[1], x[end]
(-1.0, 1.0)
```
"""
function gausslobatto(n::Integer)
    # Gauss-Legendre-Lobatto Quadrature Nodes and Weights
    if n ≤ 1
        throw(DomainError(n, "Lobatto undefined for n ≤ 1."))
    elseif n == 2
        return [-1.0, 1.0], [1.0, 1.0]
    elseif n == 3
        return [-1.0, 0.0, 1.0], [1.0 / 3, 4.0 / 3, 1.0 / 3]
    else
        # Compute via GaussJacobi:
        x, w = gaussjacobi(n - 2, 1.0, 1.0)
        @inbounds for i in 1:length(x)
            w[i] = w[i] / (1 - x[i]^2)
        end
        pushfirst!(x, -1.0)
        push!(x, 1.0)
        pushfirst!(w, 2 / (n * (n - 1)))
        push!(w, 2 / (n * (n - 1)))
        return x, w
    end
end


text:
################################################################################
### FUNCTIONS FOR EXPORTING MODEL RESULTS

function SSRecordMaker(A,N,P,D,S,T,L,Outputs,Stats; CaseDataDirectory=pwd(), PrintSpacer="*"^50)
    ################################################################################
    ### SETUP
    RecordFolder = "/records"

    # Add records folder if not present
    if !isdir(CaseDataDirectory*RecordFolder)
        mkdir(CaseDataDirectory*RecordFolder)
    end

    ################################################################################
    ### WRITE VARIABLES TO FILES

    # sl[i]
    filename = open(CaseDataDirectory*RecordFolder*"/supply.csv","w")
    print(filename, "# 1.Supply ID| 2.Node| 3.Product| 4.Variable Value")
    for i in S.ID
        print(filename, "\n"*i*"|"*S.node[i]*"|"*S.prod[i]*"|"*string(Outputs.si[i]))
    end
    close(filename)

    # dl[j]
    filename = open(CaseDataDirectory*RecordFolder*"/demand.csv","w")
    print(filename, "# 1.Demand ID| 2.Node| 3.Product| 4.Variable Value")
    for j in D.ID
        print(filename, "\n"*j*"|"*D.node[j]*"|"*D.prod[j]*"|"*string(Outputs.dj[j]))
    end
    close(filename)

    # s[n,p]
    filename = open(CaseDataDirectory*RecordFolder*"/nodal_supply.csv","w")
    print(filename, "# 1.Node| 2.Product| 3.Variable Value")
    for n in N.ID
        for p in P.ID
            print(filename, "\n"*n*"|"*p*"|"*string(Outputs.snp[n,p]))
        end
    end
    close(filename)

    # d[n,p]
    filename = open(CaseDataDirectory*RecordFolder*"/nodal_demand.csv","w")
    print(filename, "# 1.Node| 2.Product| 3.Variable Value")
    for n in N.ID
        for p in P.ID
            print(filename, "\n"*n*"|"*p*"|"*string(Outputs.dnp[n,p]))
        end
    end
    close(filename)

    # f[n,m,p]
    filename = open(CaseDataDirectory*RecordFolder*"/transport.csv","w")
    print(filename, "# 1.Arc| 2.Product| 3.Variable Value")
    for a in A.ID
        for p in P.ID
            print(filename, "\n"*a*"|"*p*"|"*string(Outputs.f[a,p]))
        end
    end
    close(filename)

    # x[t,n,q]
    filename = open(CaseDataDirectory*RecordFolder*"/consumption.csv","w")
    print(filename, "# 1.Technology| 2.Node| 3.Product| 4.Variable Value")
    for t in T.ID
        for n in N.ID
            for q in P.ID
                print(filename, "\n"*t*"|"*n*"|"*q*"|"*string(Outputs.ξcon[t,n,q]))
            end
        end
    end
    close(filename)

    # g[t,n,p]
    filename = open(CaseDataDirectory*RecordFolder*"/generation.csv","w")
    print(filename, "# 1.Technology| 2.Node| 3.Product| 4.Variable Value")
    for t in T.ID
        for n in N.ID
            for p in P.ID
                print(filename, "\n"*t*"|"*n*"|"*p*"|"*string(Outputs.ξgen[t,n,p]))
            end
        end
    end
    close(filename)

    # π[n,p]
    filename = open(CaseDataDirectory*RecordFolder*"/nodal_price.csv","w")
    print(filename, "# 1.Node| 2.Product| 3.Variable Value")
    for n in N.ID
        for p in P.ID
            print(filename, "\n"*n*"|"*p*"|"*string(Outputs.πNP[n,p]))
        end
    end
    close(filename)

    # π_transport[a,p]
    filename = open(CaseDataDirectory*RecordFolder*"/transport_price.csv","w")
    print(filename, "# 1. Arc| 2.Node 1| 3.Node 2| 4.Product| 5.Variable Value")
    for a in A.ID
        for p in P.ID
            print(filename, "\n"*a*"|"*A.node_s[a]*"|"*A.node_r[a]*"|"*p*"|"*string(Outputs.πA[a,p]))
        end
    end
    close(filename)

    # cp_tech[t,n,z]
    filename = open(CaseDataDirectory*RecordFolder*"/technology_price.csv","w")
    print(filename, "# 1.Technology| 2.Node| 3.Product| 4.Variable Value")
    for t in T.ID
        for n in N.ID
            print(filename, "\n"*t*"|"*n*"|"*string(Outputs.πT[t,n]))
        end
    end
    close(filename)

    # PhiSupply[i]
    filename = open(CaseDataDirectory*RecordFolder*"/supply_profit.csv","w")
    print(filename, "# 1.Supply ID| 2.Node| 3.Product| 4.Variable Value")
    for i in S.ID
        print(filename, "\n"*i*"|"*S.node[i]*"|"*S.prod[i]*"|"*string(Outputs.Φs[i]))
    end
    close(filename)

    # PhiDemand[i]
    filename = open(CaseDataDirectory*RecordFolder*"/demand_profit.csv","w")
    print(filename, "# 1.Demand ID| 2.Node| 3.Product| 4.Variable Value")
    for j in D.ID
        print(filename, "\n"*j*"|"*D.node[j]*"|"*D.prod[j]*"|"*string(Outputs.Φd[j]))
    end
    close(filename)

    # PhiTransport[n,m,p]
    filename = open(CaseDataDirectory*RecordFolder*"/transport_profit.csv","w")
    print(filename, "# 1. Arc| 2.Node 1| 3.Node 2| 4.Product| 5.Variable Value")
    for a in A.ID
        for p in P.ID
            print(filename, "\n"*a*"|"*A.node_s[a]*"|"*A.node_r[a]*"|"*p*"|"*string(Outputs.Φf[a,p]))
        end
    end
    close(filename)

    # PhiTech[t,n]
    filename = open(CaseDataDirectory*RecordFolder*"/technology_profit.csv","w")
    print(filename, "# 1.Technology| 2.Node| 3.Variable Value")
    for t in T.ID
        for n in N.ID
            print(filename, "\n"*t*"|"*n*"|"*string(Outputs.Φξ[t,n]))
        end
    end
    close(filename)

    # Stats
    filename = open(CaseDataDirectory*RecordFolder*"/model_stats.txt","w")
        print(filename, PrintSpacer*"\nModel statistics:")
        print(filename, "\nVariables: "*string(Stats.NumVars))
        print(filename, "\nTotal inequality constraints: "*string(Stats.TotalIneqCons))
        print(filename, "\nTotal equality constraints: "*string(Stats.TotalEqCons))
        print(filename, "\nVariable bounds: "*string(Stats.NumVarBounds))
        print(filename, "\nModel inequality constraints: "*string(Stats.ModelIneqCons))
        print(filename, "\nModel equality constraints: "*string(Stats.ModelEqCons))
        print(filename, "\n"*PrintSpacer*"\nObjective value: "*string(Outputs.z)*"\n"*PrintSpacer)
    close(filename)

    ################################################################################
    ### END OF CODE ###
    println(PrintSpacer,"\nRecords Saved!\n",PrintSpacer)
end


text:
module SpreadsheetFunctions

using DataFrames

export SUM, PRODUCT, CONCATENATE

struct Default
    value::Any
end

const NonNumber = Union{Nothing, String, Default}

export flatten_args

flatten_args(a) = Any[a]
flatten_args(a::AbstractArray) = flatten_args(a...)
function flatten_args(a, b...)
    args = flatten_args(a)
    for arg in b
        append!(args, flatten_args(arg))
    end
    return args
end

get_result(d::Default) = d.value
get_result(a) = a

function apply(f, default, a...)
    args = flatten_args(a...)
    res = foldl(f, args; init=Default(default))
    return get_result(res);
end

_SUM(a, b) = _SUM(b, a)
_SUM(a::Real, b::Real) = a + b
_SUM(a::NonNumber, b::Real) = b
_SUM(a::NonNumber, b::NonNumber) = Default(0)
_SUM(a::DataFrames.DataFrameRow, b) = SUM(Vector(a)..., b)
_SUM(a::DataFrame, b) = SUM(values(eachcol(a))..., b)

SUM(a...) = apply(_SUM, 0, a...)

_CONCATENATE(a::Union{String, Real}, b::Union{String, Real}) = string(a) * string(b)
_CONCATENATE(a::Union{String, Real}, b::Union{Nothing, Default}) = a
_CONCATENATE(a::Union{Nothing, Default}, b::Union{String, Real}) = b
_CONCATENATE(a::Union{Nothing, Default}, b::Union{Nothing, Default}) = Default("")
_CONCATENATE(a::DataFrames.DataFrameRow, b) = CONCATENATE(Vector(a)..., b)
_CONCATENATE(a::DataFrame, b) = CONCATENATE(values(eachcol(a))..., b)
_CONCATENATE(a, b::DataFrames.DataFrameRow) = CONCATENATE(b, Vector(a)...)
_CONCATENATE(a, b::DataFrame) = CONCATENATE(b, values(eachcol(a))...)

CONCATENATE(a...) = apply(_CONCATENATE, "", a...)

_PRODUCT(a, b) = _PRODUCT(b, a)
_PRODUCT(a::Real, b::Real) = a * b
_PRODUCT(a::NonNumber, b::Real) = b
_PRODUCT(a::NonNumber, b::NonNumber) = Default(0)
_PRODUCT(a::DataFrames.DataFrameRow, b) = PRODUCT(Vector(a)..., b)
_PRODUCT(a::DataFrame, b) = PRODUCT(values(eachcol(a))..., b)

PRODUCT(a...) = apply(_PRODUCT, 0, a...)

end


text:
"""
`Basics.add()`  ... 'adds' various quantities in an obvious and appropriate manner.

  + `(pota::Radial.Potential, potb::Radial.Potential)`
    ... to add two radial potentials together that are defined on the same grid. A potential::RadialPotential is returned
        that inherits its radial size from the potential that is defined in a larger range of r-values.
"""
function  Basics.add(pota::Radial.Potential, potb::Radial.Potential)
    if  pota.grid.NoPoints != potb.grid.NoPoints  ||   pota.grid.rnt != potb.grid.rnt  ||  pota.grid.h != potb.grid.h  ||
        pota.grid.hp != potb.grid.hp                   error("stop a")   end
    name = pota.name * "+" * potb.name;   nmax = max(length(pota.Zr), length(potb.Zr));   nmin = min(length(pota.Zr), length(potb.Zr))
    Zr   = zeros(nmax);
    nx = length(pota.Zr);    Zr[1:nx] = Zr[1:nx] + pota.Zr[1:nx]
    nx = length(potb.Zr);    Zr[1:nx] = Zr[1:nx] + potb.Zr[1:nx]
    Zr[nmin+1:nmax] .= Zr[nmin]

    potential = Radial.Potential(name, Zr, pota.grid)
    return( potential )
end

"""
`Basics.addZerosToCsfR(nz::Int64, csf::CsfR)`
    ... 'adds' a number of zeros to the occupation, seniority, etc. of a given CsfR.
"""
function  Basics.addZerosToCsfR(nz::Int64, csf::CsfR)
    !csf.useStandardSubshells   &&   error("Zeros can only be added to CSF with standard subshell order.")

    println("nz = $nz")
    occupation = csf.occupation;   seniority = csf.seniority;    subshellJ = csf.subshellJ;    subshellX = csf.subshellX

    for  i = 1:nz
        push!(occupation, 0);    push!(seniority, 0);    push!(subshellJ, AngularJ64(0));    push!(subshellX, subshellJ[end])
    end

    newCsf = CsfR( csf.useStandardSubshells, csf.J, csf.parity, occupation, seniority, subshellJ, subshellX, Subshell[])
    return(newCsf)
end


text:
function linear(t_0, t_n, i, n)
    return t_0 - i * (t_0 - t_n)/n
end

ln(x) = log(e, x)

function logaritmic(t_0, t_n, i, n)
    return t_0 - i ^ ( ln(t_0 - t_n) / ln(n) )
end

function hyperbolic_cos(t_0, t_n, i, n)
     return (t_0 - t_n) / cosh(6.0 * (i/n))
end

function accept_default(Δc, T, t_0)
    return exp( -Δc / T )
end


text: using Test, OffsetArrays


text:
using dio
using Test

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


text:
import Base: empty!, getindex, show
import Printf: @printf

if !isdefined(Main, :Partition)
include("core.jl")
end

"""
    entropy(partition)

Return average entropy in a partition.
The entropy of a bucket of size n is defined as lg(n), i.e. the number of additional bits
required to distinguish items in the bucket. The average is taken over all answers
in the partition.
```
julia> p=partition(["aback", "abase", "abate", "abbey"], "abbot")
ggybb aback abase
ggyby abate
gggbb abbey

julia> entropy(p)
0.5
```
"""
function entropy(p::Partition)
    sum = 0.0
    count = 0
    for i in p.nonEmptyBuckets
        n = length(p.buckets[i])
        sum += n * log(2, n)
        count += n
    end
    return sum / count
end

"""
    entropy(answers, guess)

Return remaining entropy, after partitioning possible answers with a guess.

Example:
```
julia> entropy(["aback", "abase", "abate", "abbey"], "abbot")
0.5
```
"""
entropy(answers::Vector{String}, guess::String) = entropy(partition(answers, guess))

"""
    minEntropy(partition, answers, guesses, maxCount)

Low-level variant of minEntropy. Caller passes in a Partition to use as scratch space.
Returns (entropy, guessIndex) pairs that minimize remaining entropy, in order from least to highest entropy.
"""
function minEntropy(scratch::Partition, answers::Vector{AnswerIndex}, guesses::Vector{GuessIndex}, maxCount::Integer)
    e = zeros(Float64, length(guesses))
    for i in 1:length(guesses)
        partition!(scratch, answers, guesses[i])
        e[i] = entropy(scratch)
    end
    perm = sortperm(e)
    return [(guesses[j], e[j]) for j in perm[1:min(maxCount,end)]]
end

"""
    minEntropy(answers, guesses, maxCount)

Return (entropy, guess) pairs that minimize remaining entropy, in order from least to highest entropy.
Returns min(maxCount, length(guesses)) elements.

Example:
```
julia> minEntropy(answerWords, guessWords, 4)
4-element Vector{Tuple{String, Float64}}:
 ("soare", 5.290836367768747)
 ("roate", 5.2940171538556315)
 ("raise", 5.29888678732612)
 ("raile", 5.311086768195722)
```
"""
function minEntropy(answers::Vector{String}, guesses::Vector{String}, maxCount::Integer)
    p = Partition()
    best = minEntropy(p, indicesOfWords(answers, answerWords), indicesOfWords(guesses, guessWords), maxCount)
    return [(guessWords[pair[1]], pair[2]) for pair in best]
end

"""
    search(answers, guesses, width)

Low-level overload of search that takes indices instead of strings.
"""
function search(answers::Vector{AnswerIndex}, guesses::Vector{GuessIndex}, width::Integer)
    minAvg = typemax(Float64)
    minGuess = nothing
    p = Partition()
    for (guess, entropy) in minEntropy(p, answers, guesses, width)
        partition!(p, answers, guess)
        if length(p) == 1
            # Useless guess
            continue
        end
        # sum accumulates number of guesses to guess each possible answer.
        sum = 0.0
        for i in p.nonEmptyBuckets
            # Get possible answers in the bucket
            a = p.buckets[i]
            if length(a) > 2
                sum += search(a, guesses, width)[2] * length(a)
            elseif length(a) == 2
                # Choose one of the answers as the next guess.
                # Takes one guess if lucky, two if not, for total of 3 over the two answers.
                sum += 3.0
            elseif answerWords[a[1]] == guessWords[guess]
                # Guess matched the answer
                sum += 0.0
            else
                # Only one possible answer left
                sum += 1.0
            end
        end
        avg = sum / length(answers) + 1
        if avg < minAvg
            minAvg = avg
            minGuess = guess
        end
    end
    return (minGuess, minAvg)
end

"""
    search(answers, guesses, width)

Do heuristic recursive search that minimizes average number of guesses to pick out an answer.
At each tree level, width guesses that minimize entropy are tried.

Returns pair with best first guess and average number of guesses (including that guess).
```
julia> search(["aback", "abase", "abate", "abbey"], ["aback", "abase", "abate", "abbey"], 4)
("abase", 1.75)
```
"""
function search(answers::Vector{String}, guesses::Vector{String}, width::Integer)
    (guessIndex, avg) = search(indicesOfWords(answers, answerWords), indicesOfWords(guesses, guessWords), width)
    return (guessWords[guessIndex], avg)
end

"""
    benchmark()

Time method `search` for successively wider windows.
```
julia> benchmark()
  3.843582 seconds (239.37 k allocations: 77.807 MiB, 0.23% gc time)
1 soare 3.532613
 10.000957 seconds (831.42 k allocations: 279.804 MiB, 0.33% gc time)
2 soare 3.524406
 19.578414 seconds (1.78 M allocations: 608.635 MiB, 0.31% gc time)
3 soare 3.520086
 33.192401 seconds (3.19 M allocations: 1.075 GiB, 0.32% gc time)
4 raile 3.514903
 52.062859 seconds (5.03 M allocations: 1.703 GiB, 0.33% gc time)
5 reast 3.480346
 76.379512 seconds (7.41 M allocations: 2.523 GiB, 0.33% gc time)
6 slate 3.476026
 99.450780 seconds (10.45 M allocations: 3.576 GiB, 0.33% gc time)
7 slate 3.474730
130.135821 seconds (14.07 M allocations: 4.829 GiB, 0.33% gc time)
8 slate 3.474730
```
"""
function benchmark()
    for width in 1:20
        p = @time search(answerWords, guessWords, width)
        @printf("%d %s %.6f\n", width, p[1], p[2])
    end
end

nothing

max_stars_repo_name: ArchRobison/WordleAnalysis

max_issues_repo_name: ArchRobison/WordleAnalysis

max_forks_repo_name: ArchRobison/WordleAnalysis


text:
using JLD
# using Pandas
using JSON
using Base.Profile
@everywhere using Base.Test
@everywhere import Base.convert
@everywhere using MAT
@everywhere using ForwardDiff
@everywhere using Polynomials
@everywhere import ForwardDiff.DiffBase

# local
# @everywhere include("/home/genis/cluster_home/model_fitting_julia/gradient_utils.jl")
# @everywhere include("/home/genis/cluster_home/model_fitting_julia/functions_julia.jl")
# path="/home/genis/cluster_home/model_fitting_julia/"

#portatil

path="/home/emma/fit_behaviour_julia/"
@everywhere include(path*"functions/functions_julia.jl")
@everywhere include(path*"functions/utils.jl")
@everywhere include(path*"functions/analysis.jl")
@everywhere include(path*"functions/parameter_functions.jl")

#SUBJECTS=[25,35,37]
SUBJECTS = [25]
# filename="synthetic_data_DW_Nstim10000.jld"
#model_fit="Stconst_DW"
#model_fit="StLogisticDw_2"
#model_fit="StLinTimeDw"
#model_fit="PI_Stconst"
#for model_fit in ["StLinTimePi"]
#for model_fit in ["StLinTimeDw","StLinDw","StLogDw_1"]
for model_fit in ["StLinTimeDw"]
	Nframes=20
	Tframe=0.5

	#subject=6
	for subject in SUBJECTS
		println("subject: ",subject)
		# Nconditions=1
		j=1
		path="/home/emma/fit_behaviour_julia/"
		path_files=path*"results/model_fit"*model_fit*"/subject"*string(subject)*"_1sec/"

		save_fit=analysis_minimization_subjects(path_files)
		path_save=path*"results/summary/model_fit"*model_fit*"/"
		mkpath(path_save)
		filename_save=path_save*"rat"*string(subject)*"_1sec.jld"
		JLD.save(filename_save,save_fit)
		# println(type(save_fit))
		# df=Pandas.DataFrame(save_fit)
		filename_save=path_save*"rat"*string(subject)*"_1sec.json"

		stringdata = JSON.json(save_fit)

		# write the file with the stringdata variable information
		open(filename_save, "w") do f
		        write(f, stringdata)
		     end
		# Pandas.to_json(save_fit,filename_save)
	end
end

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject