text:
#Brownian Bridge path construction
export BrownianBridgeConstruction, transform!, transformToUnit!

mutable struct BrownianBridgeConstruction
    sqrtdt::Vector{Float64}
    stdDev::Vector{Float64}
    bridgeIndex::Vector{Int}
    leftIndex::Vector{Int}
    rightIndex::Vector{Int}
    leftWeight::Vector{Float64}
    rightWeight::Vector{Float64}
    pmap::Vector{Int}
end
#=
This is a port of Quantlib code

 Copyright (C) 2003 Ferdinando Ametrano
 Copyright (C) 2006 StatPro Italia srl
 Copyright (C) 2009 Bojan Nikolic

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
=#
#===========================================================================
 The following copyright notice applies to the original code,

 Copyright (C) 2002 Peter Jäckel "Monte Carlo Methods in Finance".
 All rights reserved.

 Permission to use, copy, modify, and distribute this software is freely
 granted, provided that this notice is preserved.
 ===========================================================================#

#t[1] is first non-zero time, implicitly assumes that first time = 0
function BrownianBridgeConstruction(t::Vector{Float64})
    size = length(t)
    sqrtdt = Vector{Float64}(undef, size)
    sqrtdt[1] = sqrt(t[1])
    @inbounds for i = 2:size
        sqrtdt[i] = sqrt(t[i] - t[i-1])
    end

    # map is used to indicate which points are already constructed.
    # If map[i] is one, path point i is yet unconstructed.
    # map[i]-1 is the index of the variate that constructs
    # the path point i.
    pmap = Vector{Int}(undef, size)
    fill!(pmap, 1)
    bridgeIndex = Vector{Int}(undef, size)
    leftIndex = Vector{Int}(undef, size)
    rightIndex = Vector{Int}(undef, size)
    stdDev = Vector{Float64}(undef, size)
    leftWeight = Vector{Float64}(undef, size)
    rightWeight = Vector{Float64}(undef, size)

    #  The first point in the construction is the global step.
    pmap[size] = 2
    #  The global step is constructed from the first variate.
    bridgeIndex[1] = size
    #  The variance of the global step
    stdDev[1] = sqrt(t[size])
    #  The global step to the last point in time is special.
    leftWeight[1] = 0.0
    rightWeight[1] = 0.0
    j = 1
    for i = 2:size
        # Find the next unpopulated entry in the map.
        while pmap[j] != 1
            j += 1
        end
        k = j
        # Find the next populated entry in the map from there.
        while pmap[k] == 1
            k += 1
        end
        # l-1 is now the index of the point to be constructed next.
        l = j + ((k - 1 - j) >> 1)
        pmap[l] = i
        # The i-th Gaussian variate will be used to set point l-1.
        bridgeIndex[i] = l
        leftIndex[i] = j
        rightIndex[i] = k
        if j != 1
            leftWeight[i] = (t[k] - t[l]) / (t[k] - t[j-1])
            rightWeight[i] = (t[l] - t[j-1]) / (t[k] - t[j-1])
            stdDev[i] = sqrt(((t[l] - t[j-1]) * (t[k] - t[l])) / (t[k] - t[j-1]))
        else
            leftWeight[i] = (t[k] - t[l]) / t[k]
            rightWeight[i] = t[l] / t[k]
            stdDev[i] = sqrt(t[l] * (t[k] - t[l]) / t[k])
        end
        j = k + 1
        if j > size
            j = 1 #  wrap around
        end
    end
    return BrownianBridgeConstruction(
        sqrtdt,
        stdDev,
        bridgeIndex,
        leftIndex,
        rightIndex,
        leftWeight,
        rightWeight,
        pmap,
    )
end

function transformToUnit!(
    bb::BrownianBridgeConstruction,
    in::Vector{Float64},
    out::Vector{Float64},
)
    transform!(bb, in, out)
    size = length(bb.stdDev)
    @. out /= bb.sqrtdt
    return out
end

function transform!(
    bb::BrownianBridgeConstruction,
    in::Vector{Float64},
    out::Vector{Float64},
)
    size = length(bb.stdDev)
    # We use output to store the path...
    @inbounds out[size] = bb.stdDev[1] * in[1]
    @inbounds for i = 2:size
        j = bb.leftIndex[i]
        k = bb.rightIndex[i]
        l = bb.bridgeIndex[i]
        if j != 1
            @inbounds out[l] =
                bb.leftWeight[i] * out[j-1] +
                bb.rightWeight[i] * out[k] +
                bb.stdDev[i] * in[i]
        else
            @inbounds out[l] = bb.rightWeight[i] * out[k] + bb.stdDev[i] * in[i]
        end
    end
    # ...after which, we calculate the variations
    @inbounds for i = size:-1:2
        out[i] -= out[i-1]
    end
    return out
end

function transform!(
    bb::BrownianBridgeConstruction,
    in::Vector{Float64},  #of (d * ntimes). 1,...,d used for first time, d+1,...d+d for second time.
    out::Array{Float64,2}, #size (d,ntimes) retrieve column(ti) = out[:,ti] (column major order)
)
    d = size(out, 1) #number of rows
    bsize = length(bb.stdDev)
    # We use output to store the path...
    @inbounds for di = 1:d
        out[di, bsize] = bb.stdDev[1] * in[di]
    end
    @inbounds for i = 2:bsize
        j = bb.leftIndex[i]
        k = bb.rightIndex[i]
        l = bb.bridgeIndex[i]
        if j != 1
            for di = 1:d
                out[di, l] =
                    bb.leftWeight[i] * out[di, j-1] +
                    bb.rightWeight[i] * out[di, k] +
                    bb.stdDev[i] * in[(i-1)*d+di]
            end
        else
            for di = 1:d
                out[di, l] =
                    bb.rightWeight[i] * (out[di, k]) + bb.stdDev[i] * in[(i-1)*d+di]
            end
        end
    end
    # ...after which, we calculate the variations
    @inbounds for i = bsize:-1:2
        for di = 1:d
            out[di, i] -= out[di, i-1]
        end
    end
    return out
end


text:
push!(LOAD_PATH, "./src")
using DelimitedFiles
using Plots
using SetupDG
using UnPack
using UniformTriMesh
using PlotThemes
using ToeplitzMatrices
using Basis2DTri
using Basis1D
using LinearAlgebra

S_N(x) = @. sin(pi*x/h)/(2*pi/h)/tan(x/2)
function vandermonde_Sinc(h,r)
    N = convert(Int, 2*pi/h)
    V = zeros(length(r),N)
    for n = 1:N
        V[:,n] = S_N(r.-n*h)
    end
    V[1,1] = 1
    V[end,end] = 1
    return V
end

T = 2.0
N = 4
Nq_P = countlines("src/QuadratureData/quad_nodes_tri_N$(2*N).txt")
Nd = 5
K1D = 128
K = K1D*K1D*2
Np_F = 8
Nplot = 20
Nplot_f = 250
Mach = ".7"
delta = ".05"
extra_attribute = "varwidthsquared"

folder_name = "K$(K1D)NPF$(Np_F)N$(N)delta$(delta)Mach$(Mach)$(extra_attribute)"
file_name = "array_4"
Q = readdlm("./examples/ES_Fourier-DG/$folder_name/$file_name.txt",'\t',Float64,'\n')
xq = readdlm("./examples/ES_Fourier-DG/$folder_name/xq.txt",'\t',Float64,'\n')
yq = readdlm("./examples/ES_Fourier-DG/$folder_name/yq.txt",'\t',Float64,'\n')

xq = reshape(xq,Nq_P,Np_F,K)
yq = reshape(yq,Nq_P,Np_F,K)

h       = 2*pi/Np_F
column  = [0; .5*(-1).^(1:Np_F-1).*cot.((1:Np_F-1)*h/2)]
Dt      = Array{Float64,2}(Toeplitz(column,column[[1;Np_F:-1:2]]))

rd = init_reference_tri(N);
@unpack fv,Nfaces,r,s,VDM,V1,Vp,Dr,Ds,rf,sf,wf,nrJ,nsJ,rq,sq,wq,Vq,M,Pq,Vf,LIFT = rd
rp, sp = Basis2DTri.equi_nodes_2D(Nplot)
Vp = Basis2DTri.vandermonde_2D(N,rp,sp)/VDM
Np = length(rp)

"Mesh related variables"
# Initialize 2D triangular mesh
VX,VY,EToV = uniform_tri_mesh(K1D,K1D)
@. VX = VX
@. VY = VY
md   = init_mesh((VX,VY),EToV,rd)
# Intialize triangular prism
VX   = repeat(VX,2)
VY   = repeat(VY,2)
VZ   = [2/Np_F*ones((K1D+1)*(K1D+1),1); ones((K1D+1)*(K1D+1),1)]
@unpack rxJ,sxJ,ryJ,syJ = md

JP = 1/K1D^2
JF = 1/pi
J = JF*JP
Dr = JF*Dr
Ds = JF*Ds
Dt = JP*Dt
VpPq = Vp*Pq

Q = reshape(Q,Nq_P,Np_F,Nd,K)
u = reshape(Q[:,:,2,:],Nq_P,Np_F*K)
v = reshape(Q[:,:,3,:],Nq_P,Np_F*K)
w = reshape(Q[:,:,4,:],Nq_P,Np_F*K)

dudr = Vq*Dr*Pq*u
dvdr = Vq*Dr*Pq*v
dwdr = Vq*Dr*Pq*w
duds = Vq*Ds*Pq*u
dvds = Vq*Ds*Pq*v
dwds = Vq*Ds*Pq*w

dudt = zeros(size(dudr))
dvdt = zeros(size(dvdr))
dwdt = zeros(size(dwdr))

for k = 1:K
   f_idx = 1+(k-1)*Np_F:k*Np_F
   dudt[:,f_idx] = u[:,f_idx]*Dt'
   dvdt[:,f_idx] = v[:,f_idx]*Dt'
   dwdt[:,f_idx] = w[:,f_idx]*Dt'
end

q_criterion = zeros(Nq_P,Np_F*K)
tz = 1/JF
for k = 1:K
   rx,ry,sx,sy = rxJ[1,k],ryJ[1,k],sxJ[1,k],syJ[1,k]
   for f = 1:Np_F
      for q = 1:Nq_P
         idx = f+(k-1)*Np_F
         v_tensor = [dudr[q,idx]*rx+duds[q,idx]*sx dudr[q,idx]*ry+duds[q,idx]*sy dudt[q,idx]*tz;
                     dvdr[q,idx]*rx+dvds[q,idx]*sx dvdr[q,idx]*ry+dvds[q,idx]*sy dvdt[q,idx]*tz;
                     dwdr[q,idx]*rx+dwds[q,idx]*sx dwdr[q,idx]*ry+dwds[q,idx]*sy dwdt[q,idx]*tz;]
         s_tensor = 1/2*(v_tensor+transpose(v_tensor))
         Ω_tensor = 1/2*(v_tensor-transpose(v_tensor))
         q_criterion[q,idx] = 1/2*(norm(Ω_tensor,2)^2-norm(s_tensor,2)^2)
      end
   end
end

q_criterion = q_criterion/maximum(abs.(q_criterion))
q_criterion = reshape(q_criterion,Nq_P,Np_F,K)
zq = reshape(repeat(collect(-1+2/Np_F:(2/Np_F):1),inner=(1,Nq_P),outer=(K,1))',Nq_P,Np_F,K)
zp = LinRange(-1+2/Np_F,1,Nplot_f)
tp = LinRange(h,2*pi,Nplot_f)
VDM_F = vandermonde_Sinc(h,tp)
V2 = vandermonde_1D(1,tp)/vandermonde_1D(1,LinRange(h,2*pi,Np_F))

x_plot = zeros(size(Vp,1),Nplot_f,K)
y_plot = zeros(size(Vp,1),Nplot_f,K)
z_plot = zeros(size(Vp,1),Nplot_f,K)
q_criterion_plot = zeros(size(Vp,1),Nplot_f,K)
for k = 1:K
   x_plot[:,:,k] = Vp*Pq*xq[:,:,k]*V2'
   y_plot[:,:,k] = Vp*Pq*yq[:,:,k]*V2'
   z_plot[:,:,k] = Vp*Pq*zq[:,:,k]*V2'
   q_criterion_plot[:,:,k] = Vp*Pq*q_criterion[:,:,k]*VDM_F'
end

open("./examples/ES_Fourier-DG/q_criterion.txt","w") do io
   writedlm(io,Array(q_criterion_plot))
end

open("./examples/ES_Fourier-DG/q_criterion_x.txt","w") do io
   writedlm(io,Array(x_plot))
end

open("./examples/ES_Fourier-DG/q_criterion_y.txt","w") do io
   writedlm(io,Array(y_plot))
end

open("./examples/ES_Fourier-DG/q_criterion_z.txt","w") do io
   writedlm(io,Array(z_plot))
end

ext: jl

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
"""
The *LaplaceModule* contains all functionality to work with the Laplace transformation
1/x ≈ ∑ₜweight[t] * exp(-x node[t]).
The nodes and weights are contained within the *LaplacePoints* type. Hackbusch and coworkers
have optimized and pretabulated these for different x-ranges according to minmax
(Takatsuka, Ten-No, Hackbusch, JCP, 129 (2008), 044112).
Their pretabulated files can be obtained with *downloadLaplacePointsHackbusch*. Use *readLaplacePointsHackbusch* to
generate a LaplacePoints object from that file. After downloading, one can also use *findLaplacePointsHackbuschPretableLarger*/*findLaplacePointsHackbuschPretableSmaller* to obtain the best guesses for a given range. The ideal Laplace
problem considers x∈ [1,R]. For the transformation to the x∈ [A,B]-type Laplace problem utility functions are provided (*transformRangeToIdealLaplace*,*transformWeightFromIdealLaplace*,*transformNodeFromIdealLaplace*,*transformLaplacePointFromIdealLaplace*).
"""
module LaplaceModule
export LaplacePoints, downloadLaplacePointsHackbusch, findLaplacePointsHackbuschPretableLarger, findLaplacePointsHackbuschPretableSmaller, transformRangeToIdealLaplace, transformNodeFromIdealLaplace, transformWeightFromIdealLaplace, transformLaplacePointFromIdealLaplace
using ..DocumentationModule
using ..BaseModule
using ZipFile
using Printf
import Base.show
import Base.==

"""
LaplacePoints contain the collection of weights and corresponding nodes for the quadrature
according to
1/x = ∫ exp(-xt) ≈ ∑ₜweight[t] * exp(-x node[t])
"""
struct LaplacePoints
  weights::Array{Float64,1}
  nodes::Array{Float64,1}
end

function show(io::IO,::MIME"text/plain",lp::LaplacePoints)
  @assert length(lp.weights) == length(lp.nodes)
  println(io,"  Exponential Factor         Weight")
  println(io,"  ------------------         ------")
  for (node,weight) in zip(lp.nodes,lp.weights)
    print(io,@sprintf("  %.16f         %f\n",node,weight))
  end
end

"""
Evaluates to the name of the Laplace point file which corresponds to the given number of points and limit
in Hackbusch nomenclature. These files are contained in the zip-file obtained via downloadLaplacePointsHackbusch.
"""
function computeFilenameHackbusch(numberOfPoints::Integer,R::Real)
  k = numberOfPoints
  kstr = @sprintf("k%02d",k)
  Rstr = replace(replace(@sprintf(".%.0E",R),r"\+0*(\d+)"=>s"\1"),r"(\d+)E(\d+)"=>s"\1_\2")
  filename = "1_x$kstr$Rstr"
  return filename
end

"""
Extract the number of points (k) and the range (1,R) from the filename generated with computeFilenameHackbusch
"""
function computeFilenameHackbuschReverse(filename::AbstractString)
  m = match(r"1_xk(\d+).(\d+)_(\d+)",filename)
  k = parse(Int64,m.captures[1])
  R = parse(Float64,"$(m.captures[2])E$(m.captures[3])")
  return (k,R)
end

"""
Download a zip-file of pretabulated Laplace points to target location (default:"hackbusch")
"""
function downloadLaplacePointsHackbusch(target::AbstractString="hackbusch")
  download("https://www.mis.mpg.de/scicomp/EXP_SUM/1_x/1_xData.zip",target)
end
@doc GenericCitation("https://www.mis.mpg.de/scicomp/EXP_SUM/1_x/") downloadLaplacePointsHackbusch

"""
Generate a LaplacePoints object by extracting from the pretable zip-file obtained with downloadLaplacePointsHackbusch.
"""
function readLaplacePointsHackbusch(numberOfPoints::Integer,R::Real,zipfile::AbstractString="hackbusch")
  k = numberOfPoints
  dir = ZipFile.Reader(zipfile)
  file = ZipFile.ReadableFile
  for f in dir.files
    if f.name==computeFilenameHackbusch(k,R)
      file = f
    end
  end

  result = LaplacePoints(Array{Float64,1}(),Array{Float64,1}())
  lines = readlines(file)
  for line in 1:k
    regmatch = match(Regex(" *(?P<weight>$floatregex) .*omega.*"),lines[line])
    push!(result.weights,parse(Float64,regmatch[:weight]))
  end
  for line in k+1:2k
    regmatch = match(Regex(" *(?P<node>$floatregex) .*alpha.*"),lines[line])
    push!(result.nodes,parse(Float64,regmatch[:node]))
  end
  close(file)
  close(dir)
  return result
end

"""
Returns the LaplacePoints for the smallest range larger than the requested one for which
pretabulated Laplace points can be found. If the requested range is larger than any pretabulated one,
we expect the largest pretabulated one to be Rₖ, so the function issues a warning and
returns the LaplacePoints for that range.
"""
function findLaplacePointsHackbuschPretableLarger(numberOfPoints::Integer,R::Float64,zipfile::AbstractString="hackbusch")
  dir = ZipFile.Reader(zipfile)
  Rtight = Inf
  Rmax = 0
  for f in dir.files
    if occursin(r"1_xk(\d+).(\d+)_(\d+)",f.name)
      k,range = computeFilenameHackbuschReverse(f.name)
      if k == numberOfPoints
	Rmax = (Rmax > range) ? Rmax : range;
	Rtight = (R <= range < Rtight) ? range : Rtight;
      end
    end
  end

  if Rtight < Inf
    return (readLaplacePointsHackbusch(numberOfPoints,Rtight,zipfile),Rtight)
  else
    # Expect R > Rₖ, so just return largest available table for the given number of points
    @warn "numberOfPoints not sufficiently large ?!"
    return (readLaplacePointsHackbusch(numberOfPoints,Rmax,zipfile),Rmax)
  end
end

"""
Returns the LaplacePoints for the largest range smaller than the requested one for which
pretabulated Laplace points can be found.
"""
function findLaplacePointsHackbuschPretableSmaller(numberOfPoints::Integer,R::Float64,zipfile::AbstractString="hackbusch")
  dir = ZipFile.Reader(zipfile)
  Rtight = 0.
  for f in dir.files
    if occursin(r"1_xk(\d+).(\d+)_(\d+)",f.name)
      k,range = computeFilenameHackbuschReverse(f.name)
      if k == numberOfPoints
	Rtight = (Rtight < range <= R) ? range : Rtight;
      end
    end
  end

  if Rtight > 0
    return (readLaplacePointsHackbusch(numberOfPoints,Rtight,zipfile),Rtight)
  else
    error("No pretabulated value for $numberOfPoints points and a Range of [1., $R] or smaller.")
  end
end

"""
For the ideal Laplace problem the values x must be known to be in the range [1,R].
If 1/x can in general be within [A,B], then all nodes and weights must be scaled with 1/A,
so that x' := x/A lies within [1,B/A (=:R)].

Consider: 1/x = ∑ w exp(-α x)  =>  A/x = ∑ A w exp(-α x)
                               =>  A/x = ∑ (A w) exp(-(A α) x/A)
                   =>  1/x' = ∑ w' exp(-α' x')
"""
function transformRangeToIdealLaplace(A::Float64,B::Float64)
  return (1,B/A)
end

"""
cmp. help(transformRangeToIdealLaplace)
"""
function transformNodeFromIdealLaplace(α::Float64,A::Float64)
  return α/A
end

"""
cmp. help(transformRangeToIdealLaplace)
"""
function transformWeightFromIdealLaplace(w::Float64,A::Float64)
  return w/A
end

"""
cmp. help(transformRangeToIdealLaplace)
"""
function transformLaplacePointFromIdealLaplace(lp::LaplacePoints, A)
  LaplacePoints(map(x->transformWeightFromIdealLaplace(x,A),lp.weights), map(x->transformNodeFromIdealLaplace(x,A),lp.nodes))
end

"""
1/x ≈ ∑ᵢ wᵢ exp(-αᵢ x)
"""
function computeInverseByLaplaceApproximation(x,lp::LaplacePoints)
  sum=0.
  for (node,weight) in zip(lp.nodes,lp.weights)
    sum += weight*exp(-node*x)
  end
  return sum
end
@doc JournalCitation(["A. Takatsuka", "S. Ten-No", "W. Hackbusch"],"J. Chem. Phys.",129,044112,2008) computeInverseByLaplaceApproximation

"""
1/(x² + ω²) ≈ ∑ᵢ wᵢ sin(ω αᵢ)/ω exp(-αᵢ x)
"""
function computeRPADenominatorByDoubleLaplace(x::Float64,ω::Float64,lp::LaplacePoints)
  sum = 0.
  for (node,weight) in zip(lp.nodes,lp.weights)
    sum += weight * sin(ω*node)/ω * exp(-node*x)
  end
  return sum
end
@doc JournalCitation(["H. F. Schurkus","C. Ochsenfeld"],"J. Chem. Phys.",144,031101,2016) computeRPADenominatorByDoubleLaplace

==(x::LaplacePoints, y::LaplacePoints) = x.weights == y.weights && x.nodes == y.nodes

end # LaplaceModule


text:
##Finite Element Method Introduction

using FiniteElementDiffEq,  DiffEqPDEBase

### Setup
dx = 1//2^(3)
mesh = notime_squaremesh([0 1 0 1],dx,:dirichlet)

f = (x) -> sin.(2π.*x[:,1]).*cos.(2π.*x[:,2])
prob = PoissonProblem(f,mesh)

sol = solve(prob,FEMDiffEqPoisson())

sol = solve(prob,FEMDiffEqPoisson(),solver=:CG)

sol = solve(prob,FEMDiffEqPoisson(),solver=:GMRES)

#plot(sol)

### Test Results

#The test will only pass if the calculated L2 error is below 1e-4.
true


text:
# tid = parse(Int64, ARGS[1])
force_scale = 6.0 #50
tid = 203
fiber_size = 2
porder = 2
if length(ARGS) == 4
   global tid = parse(Int64, ARGS[1])
   global force_scale = parse(Float64, ARGS[2])
   global fiber_size = parse(Int64, ARGS[3])
   global porder = parse(Int64, ARGS[4])
end
printstyled("force_scale=$force_scale, tid=$tid, fiber_size=$fiber_size, porder=$porder\n", color=:green)
include("CommonFuncs.jl")
np = pyimport("numpy")
"""
Property:
The matrix is made of titanium, which are assumed to be elasto-plastic titanium material,
time ms
ρ = 4.5 g/cm^3;  E = 100GPa = 	100*10^9 Kg/m/s^2 = 10^6 g/cm/ms^2, K=10e+9Pa =  10^5 g/cm/ms^2
ν =0.2   σY=970 MPa = 9700

Fibers are made of SiC, which are assumed to be isotropic and elastic, with
https://www.azom.com/properties.aspx?ArticleID=42
ρ = 3.2 g/cm^3  E = 400GPa =  4*10^6 g/cm/ms^2  ν = 0.35
length scale cm
"""
prop1 = Dict("name"=> "PlaneStressPlasticity","rho"=> 4.5, "E"=> 1e+6, "nu"=> 0.2,
"sigmaY"=>0.97e+4, "K"=>1e+5)
prop2 = Dict("name"=> "PlaneStress", "rho"=> 3.2, "E"=>4e6, "nu"=>0.35)

# prop1 = prop2
ps1 = PlaneStress(prop1); H1 = ps1.H
ps2 = PlaneStress(prop2); H2 = ps2.H

T = 200.0
NT = 200

nxc, nyc = 80,40
nx, ny =  nxc*fiber_size, nyc*fiber_size
#Type 1=> SiC(fiber), type 0=>Ti(matrix), each fiber has size is k by k
fiber_fraction = 0.25
fiber_distribution = "Uniform"
ele_type = generateEleType(nxc, nyc, fiber_size, fiber_fraction, fiber_distribution)
@show "fiber fraction: ",  sum(ele_type)/(nx*ny)
# matshow(ele_type)
# savefig("test.png")
# error()

nodes, EBC, g, gt, FBC, fext, ft, npoints, node_to_point = BoundaryCondition(tid, nx, ny, porder; force_scale=force_scale)
elements = []
for j = 1:ny
    for i = 1:nx
        n = (nx*porder+1)*(j-1)*porder + (i-1)porder+1
        #element (i,j)
        if porder == 1
            #   4 ---- 3
            #
            #   1 ---- 2

            elnodes = [n, n + 1, n + 1 + (nx + 1), n + (nx + 1)]
        elseif porder == 2
            #   4 --7-- 3
            #   8   9   6
            #   1 --5-- 2
            elnodes = [n, n + 2, n + 2 + 2*(2*nx+1),  n + 2*(2*nx+1), n+1, n + 2 + (2*nx+1), n + 1 + 2*(2*nx+1), n + (2*nx+1), n+1+(2*nx+1)]
        else
            error("polynomial order error, porder= ", porder)
        end

        coords = nodes[elnodes,:]
        # 0=> matrix, 1=> fiber
        prop = ele_type[i,j] == 0 ? prop1 : prop2
        push!(elements,SmallStrainContinuum(coords,elnodes, prop, 3))
    end
end

ndofs = 2
domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)
setGeometryPoints!(domain, npoints, node_to_point)

state = zeros(domain.neqs)
∂u = zeros(domain.neqs)
globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)

assembleMassMatrix!(globdat, domain)
updateStates!(domain, globdat)

Δt = T/NT

ρ_oo = 0.0
αm = (2*ρ_oo - 1)/(ρ_oo + 1)
αf = ρ_oo/(ρ_oo + 1)

for i = 1:NT
    @info i, "/" , NT

    solver = NewmarkSolver(Δt, globdat, domain, αm, αf, 1e-4, 1e-6, 10)
    #solver = NewmarkSolver(Δt, globdat, domain, 0.5, 0.5, 1e-4, 1e-6, 10)
    # close("all")
    # visσ(domain)
    # axis("equal")
    # savefig("Debug/$i.png")
    # error()
    # close("all")
    # visσ(domain)
    # # visσ(domain,-1.5e9, 4.5e9)
    # savefig("Debug/terminal$(tid)i=$i.png")
    # if i==75
    #     close("all")
    #     visσ(domain)
    #     # visσ(domain,-1.5e9, 4.5e9)
    #     savefig("Debug/terminal$(tid)i=75.png")
    # end
end

# error()
# todo write data
if !isdir("$(@__DIR__)/Data/order$porder")
    mkdir("$(@__DIR__)/Data/order$porder")
end
if !isdir("$(@__DIR__)/Debug/order$porder")
    mkdir("$(@__DIR__)/Debug/order$porder")
end

write_data("$(@__DIR__)/Data/order$porder/$(tid)_$(force_scale)_$(fiber_size).dat", domain)
@save "Data/order$porder/domain$(tid)_$(force_scale)_$(fiber_size).jld2" domain

close("all")
visσ(domain)
axis("equal")
savefig("Debug/order$porder/terminal$(tid)_$(force_scale)_$(fiber_size).png")

close("all")
ux = [reshape(domain.history["state"][i][1:(nx*porder+1)*(ny*porder+1)], ny*porder+1, nx*porder+1)[1,end] for i = 1:length(domain.history["state"])]
plot(ux)
savefig("Debug/order$porder/ux$(tid)_$(force_scale)_$(fiber_size).png")

close("all")
uy = [reshape(domain.history["state"][i][(nx*porder+1)*(ny*porder+1)+1:end], ny*porder+1, nx*porder+1)[1,end] for i = 1:length(domain.history["state"])]
plot(uy)
savefig("Debug/order$porder/uy$(tid)_$(force_scale)_$(fiber_size).png")


text:
using ImageSegmentation, Images, Test, SimpleWeightedGraphs, LightGraphs, StaticArrays, DataStructures
using RegionTrees: isleaf, Cell, split!
using MetaGraphs: MetaGraph, clear_props!, get_prop, has_prop, set_prop!, props, vertices

using Documenter
doctest(ImageSegmentation, manual = false)

include("core.jl")
include("region_growing.jl")
include("felzenszwalb.jl")
include("fast_scanning.jl")
include("watershed.jl")
include("region_merging.jl")
include("meanshift.jl")
include("merge_segments.jl")


text:
module Parser
using ..Types
export autoparse, configparser, argparser!, shareprms!

function str2list(s, splitter=";")
    ret = []
    for e in split(s, splitter)
        if e!=""
            append!(ret, [autoparse(e)])
        end
    end
    return ret
end

# auto conversion of string
function autoparse(s)
    if s==""
        return nothing
    elseif contains(s, ";")
        return str2list(s, ";")
    elseif length(s)>1 && s[1]=='('
        @assert s[end]==')'
        # a tuple
        return tuple(autoparse(s[2:end-1])...)
    elseif contains(s, ",")
        return str2list(s, ",")
    elseif s=="yes" || s=="Yes"|| s=="y" || s=="Y" || s=="true" || s=="True"
        return true
    elseif s=="no" || s=="No" || s=="n" || s=="N" || s=="false" || s=="False"
        return false
    elseif contains(s, "/") || occursin(r"^[A-z]", s) || typeof(Meta.parse(s))==Symbol || typeof(Meta.parse(s)) == Expr
        # directory or containing alphabet not all number
        return s
    else
        return Meta.parse(s)
    end
end

"""
parse the lines
`Inputs:`
lines: cofiguration file name or lines of configuration file

`Outputs:`
pd: Dict, dictionary containing parameters
"""
function configparser(fname::AbstractString)
    if contains(fname, "[") && contains(fname, "]") && contains(fname, "=")
        # this is actually the content of the config file!
        lines = split(fname, "\n")
    else
        # this is really a file name
        lines = readlines(fname)
    end
    return configparser(lines)
end

function configparser(lines::Vector)
    # initialize the parameter dictionary
    pd = ParamDict()
    # default section name
    sec = :section
    # analysis the lines
    for l in lines
        # remove space and \n
        l = replace(l, "\n", "")
        l = replace(l, " ", "")
        if occursin(r"^\s*#", l) || occursin(r"^\s*\n", l)
            continue
        elseif occursin(r"^\s*\[.*\]", l)
            # update the section name
            m = match(r"\[.*\]", l)
            sec = Symbol( m.match[2:end-1] )
            pd[sec] = Dict()
        elseif occursin(r"^\s*.*\s*=", l)
            k, v = split(l, '=')
            k = Symbol(k)
            # assign value to dictionary
            pd[sec][k] = autoparse( v )
        end
    end
    return pd
end

"""
parse configuration file for Directed Acyclic Graph (DAG) computation
"""
function dagparser(fname::AbstractString)
    lines = readlines(fname)
    return dagparser(lines)
end

function dagparser(lines::Vector)
    # initialize the parameter dictionary
    pd = Vector(Dict)
    # temporary dictionary

    # default section name
    sec = :section
    # analysis the lines
    for l in lines
        # remove space and \n
        l = replace(l, "\n", "")
        l = replace(l, " ", "")
        if occursin(r"^\s*#", l) || occursin(r"^\s*\n", l)
            continue
        elseif occursin(r"^\s*\[.*\]", l)
            # update the section name
            m = match(r"\[.*\]", l)
            sec = Symbol( m.match[2:end-1] )
            pd[sec] = Dict()
        elseif occursin(r"^\s*.*\s*=", l)
            k, v = split(l, '=')
            k = Symbol(k)
            # assign value to dictionary
            pd[sec][k] = autoparse( v )
        end
    end
    return pd
end

"""
parameter dictionary
input can be some default dictionary

Note that all the key was one character, which is the first non '-' character!
suppose that the argument is as follows:
--flbl label.h5 --range 2,4-6
the returned dictionary will be
pd[:f] = "label.h5"
pd[:r] = [2,4,5,6]
"""
function argparser!(pd::Dict{Symbol, Any}=Dict() )
    println("default parameters: $(pd)")
    # argument table, two rows
    @assert length(ARGS) % 2 ==0
    argtbl = reshape(ARGS, 2,Int64( length(ARGS)/2))

    # traverse all the argument table columns
    for c in 1:size(argtbl,2)
        @assert argtbl[1,c][1]=='-'
        # key and value
        k = Symbol( replace(argtbl[1,c],"-","")[1] )
        v = autoparse( argtbl[2,c] )
        pd[k] = v
    end
    println("parameters after command line parsing: $(pd)")
end

# share the gneral parameters in each section
function shareprms!(pd::ParamDict, gnkey::Symbol=:gn, is_keep=true)
    @assert haskey(pd, gnkey)
    for k1 in keys(pd)
        if k1 != gnkey
            for (k2,v2) in pd[gnkey]
                pd[k1][k2] = v2
            end
        end
    end
    if !is_keep
        delete!(pd, gnkey)
    end
    return pd
end

end # end of module


text:
##################### Question 1 #####################

function bisect_left(A, p, r, v)
    i = p
    if p < r
       q = floor(Int, (p+r)/2)
       if v <= A[q]
           i = bisect_left(A, p, q, v)
       else
           i = bisect_left(A, q+1, r, v)
       end
    end
    return i
end

function bisect_right(A, p, r, v)
    i = p
    if p < r
       q = floor(Int, (p+r)/2)
       if v < A[q]
           i = bisect_right(A, p, q, v)
       else
           i = bisect_right(A, q+1, r, v)
       end
    end
    return i
end

##################### Question 2 #####################

function bisect(A, p, r, v)
    while (p <= r)
        q = floor(Int, (p+r/2))
        if (v == A[q]) return q
        elseif (v < A[q]) r = q-1
        else p = q+1
        end
    end
    return nothing
end

function qsort!(A, lo, hi) # lo - low, hi - high
    i, j = lo, hi
    while i < hi # iterate from "low" to "high" in relevant partition
        pivot = A[(lo+hi)>>>1] # logical shift (low+high) bits right to divide them by 2
        while i <= j
            while A[i] < pivot; i = i+1; end # find index less than pivot
            while A[j] > pivot; j = j-1; end # find index bigger than pivot
            if i <= j
                A[i], A[j] = A[j], A[i] # change places of the lesser and bigger element
                i, j = i+1, j-1 # next index towards middle
            end
        end
        if lo < j; qsort!(A,lo,j); end # call qsort on bottom partition
        lo, j = i, hi # start over again with the top partition
    end
    return A
end

function algdat_sort!(A) # quicksort
    p = 1
    r = length(A)
    qsort!(A, p, r)
end

function find_median(A, lower, upper)
    lo = bisect_left(A, 1, length(A)+1, lower) # lower index of interval
    hi = bisect_right(A, 1, length(A)+1, upper) - 1 # upper index of interval (-1 to return the actual last idx containing this value, not to the next)

    ctr = (lo+hi)/2 # float or integer if size of array is even or odd
    ctrLeft = (floor(Int, ctr))

    if (ctrLeft == ctr) # if array is odd, ctr index is an integer and its value the median (comparing floats)
        return A[Int(ctr)]
    else # else find median of the values to the left and right of the floating idx value
        return (A[ctrLeft] + A[ctrLeft+1]) / 2
    end
end

max_stars_repo_name: evenlwanvik-student/algdat

max_issues_repo_name: evenlwanvik-student/algdat

max_forks_repo_name: evenlwanvik-student/algdat


text:
@testset "solve_btpde.jl" begin

end


text:
module RestrictedBoltzmannMachines

import Flux
import Zygote
import ChainRulesCore

using Base: front, tail
using Random: AbstractRNG, GLOBAL_RNG, randexp, randn!, rand!, shuffle!
using LinearAlgebra: Diagonal, logdet, I, dot, norm
using Statistics: mean, std
using EllipsisNotation: (..)
using LogExpFunctions: softmax, logsumexp, log1pexp, logistic, logaddexp, logsubexp
using SpecialFunctions: erf, erfcx, logerfcx
using FillArrays: Fill, Zeros, Ones, Trues, Falses
using Flux: ADAM, Descent

include("util/util.jl")
include("util/onehot.jl")
include("util/linalg.jl")
include("util/truncated_normal.jl")
include("util/minibatches.jl")

include("layers/abstractlayer.jl")
include("layers/binary.jl")
include("layers/spin.jl")
include("layers/potts.jl")
include("layers/gaussian.jl")
include("layers/relu.jl")
include("layers/drelu.jl")
include("layers/prelu.jl")
include("layers/xrelu.jl")
include("layers/common.jl")

include("rbm.jl")

include("rbms/hopfield.jl")
include("rbms/binary.jl")
include("rbms/gaussian.jl")

include("gauge/zerosum.jl")
include("gauge/rescale_hidden.jl")

include("pseudolikelihood.jl")
include("partition.jl")
include("ais.jl")

include("centered_gradient.jl")

include("train/initialization.jl")
include("train/cd.jl")
include("train/cdad.jl")
include("train/pcd.jl")
include("train/fast.jl")
include("train/rdm.jl")
include("train/optim.jl")

include("regularize.jl")

include("metropolis.jl")

end # module


text:
const GQE{C,V} = GenericQuadExpr{C,V}
const _VariableQuadExpr{C} = GenericQuadExpr{C, VariableRef}
const _DecisionQuadExpr{C} = GenericQuadExpr{C, DecisionRef}
const _KnownQuadExpr{C} = GenericQuadExpr{C, KnownRef}
const _VQE = _VariableQuadExpr{Float64}
const _DQE = _DecisionQuadExpr{Float64}
const _KQE = _KnownQuadExpr{Float64}

struct DecisionCrossTerm
    decision::DecisionRef
    variable::VariableRef
end
Base.first(t::DecisionCrossTerm) = t.decision
Base.last(t::DecisionCrossTerm) = t.variable
Base.hash(t::DecisionCrossTerm, h::UInt) = hash(hash(t.decision) + hash(t.variable), h)
function Base.isequal(t1::DecisionCrossTerm, t2::DecisionCrossTerm)
    return (t1.decision == t2.decision && t1.variable == t2.variable)
end

struct KnownCrossTerm{V}
    known::KnownRef
    variable::V
end
Base.first(t::KnownCrossTerm) = t.known
Base.last(t::KnownCrossTerm) = t.variable
KnownVariableCrossTerm = KnownCrossTerm{VariableRef}
KnownDecisionCrossTerm = KnownCrossTerm{DecisionRef}

Base.hash(t::KnownCrossTerm, h::UInt) = hash(hash(t.known) + hash(t.variable), h)
function Base.isequal(t1::KnownCrossTerm, t2::KnownCrossTerm)
    return (t1.known == t2.known && t1.variable == t2.variable)
end

mutable struct DecisionQuadExpr{C} <: JuMP.AbstractJuMPScalar
    variables::GenericQuadExpr{C, VariableRef}
    decisions::GenericQuadExpr{C, DecisionRef}
    knowns::GenericQuadExpr{C, KnownRef}
    cross_terms::OrderedDict{DecisionCrossTerm, C}
    known_variable_terms::OrderedDict{KnownVariableCrossTerm, C}
    known_decision_terms::OrderedDict{KnownDecisionCrossTerm, C}
end
const DQE = DecisionQuadExpr{Float64}

DQE() = zero(DQE{Float64})

mutable struct BilinearDecisionExpr{C} <: JuMP.AbstractJuMPScalar
    variables::GenericAffExpr{C, VariableRef}
    decisions::GenericAffExpr{C, DecisionRef}
    knowns::GenericQuadExpr{C, KnownRef}
    known_variable_terms::OrderedDict{KnownVariableCrossTerm, C}
    known_decision_terms::OrderedDict{KnownDecisionCrossTerm, C}
end
const BDE = BilinearDecisionExpr{Float64}

# Base overrides #
# ========================== #
function Base.iszero(quad::DecisionQuadExpr)
    return iszero(quad.variables) &&
        iszero(quad.decisions) &&
        iszero(quad.knowns)
end
function Base.zero(::Type{DecisionQuadExpr{C}}) where C
    return DecisionQuadExpr{C}(
        zero(_VariableQuadExpr{C}),
        zero(_DecisionQuadExpr{C}),
        zero(_KnownQuadExpr{C}),
        OrderedDict{DecisionCrossTerm, C}(),
        OrderedDict{KnownVariableCrossTerm, C}(),
        OrderedDict{KnownDecisionCrossTerm, C}())
end
function Base.one(::Type{DecisionQuadExpr{C}}) where C
    return DecisionQuadExpr{C}(
        one(_VariableQuadExpr{C}),
        zero(_DecisionQuadExpr{C}),
        zero(_KnownQuadExpr{C}),
        OrderedDict{DecisionCrossTerm, C}(),
        OrderedDict{KnownVariableCrossTerm, C}(),
        OrderedDict{KnownDecisionCrossTerm, C}())
end
Base.zero(quad::DecisionQuadExpr) = zero(typeof(quad))
Base.one(quad::DecisionQuadExpr) =  one(typeof(quad))
function Base.copy(quad::DecisionQuadExpr{C}) where C
    return DecisionQuadExpr{C}(
        copy(quad.variables),
        copy(quad.decisions),
        copy(quad.knowns),
        copy(quad.cross_terms),
        copy(quad.known_variable_terms),
        copy(quad.known_decision_terms))
end
Base.broadcastable(quad::DecisionQuadExpr) = Ref(quad)

function Base.convert(::Type{DecisionQuadExpr{C}}, c::Number) where C
    return convert(DecisionQuadExpr{C}, convert(DecisionAffExpr{C}, c))
end
function Base.convert(::Type{DecisionQuadExpr{C}}, expr::Union{VariableRef, _VariableAffExpr}) where C
    return convert(DecisionQuadExpr{C}, convert(DecisionAffExpr{C}, expr))
end
function Base.convert(::Type{DecisionQuadExpr{C}}, expr::_VariableQuadExpr) where C
    return DecisionQuadExpr{C}(
        expr,
        zero(_DecisionQuadExpr{C}),
        zero(_KnownQuadExpr{C}),
        OrderedDict{DecisionCrossTerm, C}(),
        OrderedDict{KnownVariableCrossTerm, C}(),
        OrderedDict{KnownDecisionCrossTerm, C}())
end
function Base.convert(::Type{DecisionQuadExpr{C}}, expr::Union{DecisionRef, _DecisionAffExpr}) where C
    return convert(DecisionQuadExpr{C}, convert(DecisionAffExpr{C}, expr))
end
function Base.convert(::Type{DecisionQuadExpr{C}}, expr::_DecisionQuadExpr) where C
    return DecisionQuadExpr{C}(
        zero(_VariableQuadExpr{C}),
        expr,
        zero(_KnownQuadExpr{C}),
        OrderedDict{DecisionCrossTerm, C}(),
        OrderedDict{KnownVariableCrossTerm, C}(),
        OrderedDict{KnownDecisionCrossTerm, C}())
end
function Base.convert(::Type{DecisionQuadExpr{C}}, expr::Union{KnownRef, _KnownAffExpr}) where C
    return convert(DecisionQuadExpr{C}, convert(DecisionAffExpr{C}, expr))
end
function Base.convert(::Type{DecisionQuadExpr{C}}, expr::_KnownQuadExpr) where C
    return DecisionQuadExpr{C}(
        zero(_VariableQuadExpr{C}),
        zero(_DecisionQuadExpr{C}),
        expr,
        OrderedDict{DecisionCrossTerm, C}(),
        OrderedDict{KnownVariableCrossTerm, C}(),
        OrderedDict{KnownDecisionCrossTerm, C}())
end
function Base.convert(::Type{DecisionQuadExpr{C}}, aff::DecisionAffExpr{C}) where C
    return DecisionQuadExpr{C}(
        convert(_VariableQuadExpr{C}, aff.variables),
        convert(_DecisionQuadExpr{C}, aff.decisions),
        convert(_KnownQuadExpr{C}, aff.knowns),
        OrderedDict{DecisionCrossTerm, C}(),
        OrderedDict{KnownVariableCrossTerm, C}(),
        OrderedDict{KnownDecisionCrossTerm, C}())
end

function Base.isequal(quad::DecisionQuadExpr{C}, other::DecisionQuadExpr{C}) where C
    return isequal(quad.variables, other.variables) &&
        isequal(quad.decisions, other.decisions) &&
        isequal(quad.knowns, other.knowns) &&
        isequal(quad.cross_terms, other.cross_terms) &&
        isequal(quad.known_variable_terms, other.known_variable_terms) &&
        isequal(quad.known_decision_terms, other.known_decision_terms)
end

Base.hash(quad::DecisionQuadExpr, h::UInt) =
    hash(hash(quad.variables, h),
         hash(quad.decisions, h),
         hash(quad.knowns, h),
         hash(quad.cross_terms, h),
         hash(quad.known_variable_terms, h),
         hash(quad.known_decision_terms, h))

function sizehint!(quad::DecisionQuadExpr, n::Int, ::Type{_VAE})
    Base.sizehint!(quad.variables.aff, n)
end
function sizehint!(quad::DecisionQuadExpr, n::Int, ::Type{_DAE})
    Base.sizehint!(quad.decisions.aff, n)
end
function sizehint!(quad::DecisionQuadExpr, n::Int, ::Type{_KAE})
    Base.sizehint!(quad.knowns.aff, n)
end

# JuMP overrides #
# ========================== #
function JuMP.drop_zeros!(quad::DecisionQuadExpr)
    JuMP.drop_zeros!(quad.variables)
    JuMP.drop_zeros!(quad.decisions)
    JuMP.drop_zeros!(quad.knowns)
    dropzeros!(terms) = begin
        for (key, coef) in terms
            if iszero(coef)
                delete!(terms, key)
            end
        end
    end
    dropzeros!(quad.cross_terms)
    dropzeros!(quad.known_variable_terms)
    dropzeros!(quad.known_decision_terms)
    return nothing
end

function JuMP._affine_coefficient(f::DecisionQuadExpr, variable::VariableRef)
    return JuMP._affine_coefficient(f.variables, variable)
end
function JuMP._affine_coefficient(f::DecisionQuadExpr, decision::DecisionRef)
    return JuMP._affine_coefficient(f.decisions, decision)
end
function JuMP._affine_coefficient(f::DecisionQuadExpr, known::KnownRef)
    return JuMP._affine_coefficient(f.knowns, known)
end

function JuMP.map_coefficients_inplace!(f::Function, quad::DecisionQuadExpr)
    JuMP.map_coefficients_inplace!(f, quad.variables)
    JuMP.map_coefficients_inplace!(f, quad.decisions)
    JuMP.map_coefficients_inplace!(f, quad.knowns)
    _map_cross_terms!(quad.cross_terms)
    _map_cross_terms!(quad.known_variable_terms)
    _map_cross_terms!(quad.known_decision_terms)
    return quad
end

function _map_cross_terms!(f::Function, terms)
    for (key, value) in terms
        terms[key] = f(value)
    end
    return nothing
end

function _map_cross_terms(f::Function, terms)
    res = copy(terms)
    _map_cross_terms!(f, res)
    return res
end

function JuMP.map_coefficients(f::Function, quad::DecisionQuadExpr)
    return JuMP.map_coefficients_inplace!(f, copy(quad))
end

function JuMP.value(quad::DecisionQuadExpr, value_func::Function)
    varvalue = JuMP.value(quad.variables, value_func)
    dvarvalue = JuMP.value(quad.decisions, value_func)
    kvarvalue = JuMP.value(quad.knowns, value_func)
    T = promote_type(typeof(varvalue), typeof(dvarvalue), typeof(kvarvalue))
    ret = convert(T, varvalue) + convert(T, dvarvalue) + convert(T, kvarvalue)
    val(terms, T) = begin
        ret = zero(T)
        for (vars, coef) in quad.cross_terms
            ret += coef * value_func(vars.a) * value_func(vars.b)
        end
        ret
    end
    ret += val(quad.cross_terms, T)
    ret += val(quad.known_variable_terms, T)
    ret += val(quad.known_decision_terms, T)
    return ret
end

function JuMP.constant(quad::DecisionQuadExpr)
    return quad.variables.constant
end

function JuMP.linear_terms(quad::DecisionQuadExpr, ::Type{_VAE})
    JuMP.linear_terms(quad.variables)
end
function JuMP.linear_terms(quad::DecisionQuadExpr, ::Type{_DAE})
    JuMP.linear_terms(quad.decisions)
end
function JuMP.linear_terms(quad::DecisionQuadExpr, ::Type{_KAE})
    JuMP.linear_terms(quad.knowns)
end

function SparseArrays.dropzeros(quad::DecisionQuadExpr{C}) where C
    variables = SparseArrays.dropzeros(quad.variables)
    decisions = SparseArrays.dropzeros(quad.decisions)
    knowns = SparseArrays.dropzeros(quad.knowns)
    dropzeros!(terms) = begin
        for (key, value) in terms
            if iszero(value)
                delete!(terms, key)
            end
        end
    end
    cross_terms = copy(quad.cross_terms)
    dropzeros!(cross_terms)
    known_variable_terms = copy(quad.known_variable_terms)
    dropzeros!(known_variable_terms)
    known_decision_terms = copy(quad.known_decision_terms)
    dropzeros!(known_decision_terms)
    return DecisionQuadExpr(variables, decisions, knowns, cross_terms, known_variable_terms, known_decision_terms)
end

function JuMP._assert_isfinite(quad::Union{_DecisionQuadExpr, _KnownQuadExpr})
    JuMP._assert_isfinite(quad.aff)
    for (coef, var1, var2) in quad_terms(quad)
        isfinite(coef) || error("Invalid coefficient $coef on quadratic term $var1*$var2.")
    end
    return nothing
end

function JuMP._assert_isfinite(quad::DecisionQuadExpr)
    JuMP._assert_isfinite(quad.variables)
    JuMP._assert_isfinite(quad.decisions)
    JuMP._assert_isfinite(quad.knowns)
    _assert_isfinite(terms) = begin
        for (vars, coef) in terms
            isfinite(coef) || error("Invalid coefficient $coef on quadratic cross term $first(vars)*$last(vars).")
        end
    end
    _assert_isfinite(quad.cross_terms)
    _assert_isfinite(quad.known_variable_terms)
    _assert_isfinite(quad.known_decision_terms)
    return nothing
end

function JuMP.check_belongs_to_model(quad::DecisionQuadExpr, model::AbstractModel)
    JuMP.check_belongs_to_model(quad.variables, model)
    JuMP.check_belongs_to_model(quad.decisions, model)
    JuMP.check_belongs_to_model(quad.knowns, model)
    check_belongs_to_model(terms) = begin
        for variable_pair in keys(terms)
            JuMP.check_belongs_to_model(first(variable_pair), model)
            JuMP.check_belongs_to_model(last(variable_pair), model)
        end
    end
    check_belongs_to_model(quad.cross_terms)
    check_belongs_to_model(quad.known_variable_terms)
    check_belongs_to_model(quad.known_decision_terms)
    return nothing
end

function JuMP.function_string(mode, quad::DecisionQuadExpr, show_constant=true)
    # Decisions
    ret = ""
    decision_terms = JuMP.function_string(mode, quad.decisions)
    if decision_terms != "0"
        ret = decision_terms
    end
    known_decision_terms = _cross_terms_function_string(mode, quad.known_decision_terms)
    ret = _add_terms(ret, known_decision_terms)
    cross_terms = _cross_terms_function_string(mode, quad.cross_terms)
    ret = _add_terms(ret, cross_terms)
    # Knowns
    known_terms = JuMP.function_string(mode, quad.knowns)
    ret = _add_terms(ret, known_terms)
    known_variable_terms = _cross_terms_function_string(mode, quad.known_variable_terms)
    ret = _add_terms(ret, known_variable_terms)
    # Variables
    variable_terms = JuMP.function_string(mode, quad.variables)
    if variable_terms != "0"
        if ret == ""
            ret = variable_terms
        else
            if variable_terms[1] == '-'
                ret = ret * " - " * variable_terms[2:end]
            else
                ret = ret * " + " * variable_terms
            end
        end
    else
        if ret == ""
            return variable_terms
        end
    end
    return ret
end

function _cross_terms_function_string(mode, terms)
    term_str = Array{String}(undef, 2 * length(terms))
    elm = 1
    if length(term_str) > 0
        for (vars, coef) in terms
            pre = JuMP._is_one_for_printing(coef) ? "" : JuMP._string_round(abs(coef)) * " "

            x = function_string(mode, first(vars))
            y = function_string(mode, last(vars))

            term_str[2 * elm - 1] = JuMP._sign_string(coef)
            term_str[2 * elm] = "$pre$x"
            term_str[2 * elm] *= string(JuMP._math_symbol(mode, :times), y)
            if elm == 1
                # Correction for first term as there is no space
                # between - and variable coefficient/name
                term_str[1] = coef < zero(coef) ? "-" : ""
            end
            elm += 1
        end
    end
    ret = join(term_str[1 : 2 * (elm - 1)])
    return ret
end

# With one factor.
function JuMP.add_to_expression!(quad::DQE, other::Number)
    JuMP.add_to_expression!(quad.variables, other)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_var::VariableRef)
    JuMP.add_to_expression!(quad.variables, new_var)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_dvar::DecisionRef)
    JuMP.add_to_expression!(quad.decisions, new_dvar)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_kvar::KnownRef)
    JuMP.add_to_expression!(quad.knowns, new_kvar)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_vae::_VAE)
    JuMP.add_to_expression!(quad.variables, new_vae)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_dae::_DAE)
    JuMP.add_to_expression!(quad.decisions, new_dae)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_kae::_KAE)
    JuMP.add_to_expression!(quad.knowns, new_kae)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_vqe::_VQE)
    JuMP.add_to_expression!(quad.variables, new_vqe)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_dqe::_DQE)
    JuMP.add_to_expression!(quad.decisions, new_dqe)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_kqe::_KQE)
    JuMP.add_to_expression!(quad.knowns, new_kqe)
    return quad
end

function JuMP.add_to_expression!(lhs_quad::DQE, rhs_aff::DAE)
    JuMP.add_to_expression!(lhs_quad.variables, rhs_aff.variables)
    JuMP.add_to_expression!(lhs_quad.decisions, rhs_aff.decisions)
    JuMP.add_to_expression!(lhs_quad.knowns, rhs_aff.knowns)
    return lhs_quad
end

function JuMP.add_to_expression!(lhs_quad::DQE, rhs_quad::DQE)
    JuMP.add_to_expression!(lhs_quad.variables, rhs_quad.variables)
    JuMP.add_to_expression!(lhs_quad.decisions, rhs_quad.decisions)
    JuMP.add_to_expression!(lhs_quad.knowns, rhs_quad.knowns)
    add_cross_terms!(terms, other) = begin
        for (key, term_coef) in other
            JuMP._add_or_set!(terms, key, term_coef)
        end
    end
    add_cross_terms!(lhs_quad.cross_terms, rhs_quad.cross_terms)
    add_cross_terms!(lhs_quad.known_variable_terms, rhs_quad.known_variable_terms)
    add_cross_terms!(lhs_quad.known_decision_terms, rhs_quad.known_decision_terms)
    return lhs_quad
end

# With two factors.
function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant, new_var::VariableRef)
    JuMP.add_to_expression!(quad.variables, new_coef, new_var)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_var::VariableRef, new_coef::_Constant)
    JuMP.add_to_expression!(quad.variables, new_coef, new_var)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant, new_dvar::DecisionRef)
    JuMP.add_to_expression!(quad.decisions, new_coef, new_dvar)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_dvar::DecisionRef, new_coef::_Constant)
    JuMP.add_to_expression!(quad.decisions, new_coef, new_dvar)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant, new_kvar::KnownRef)
    JuMP.add_to_expression!(quad.knowns, new_coef, new_kvar)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_kvar::KnownRef, new_coef::_Constant)
    JuMP.add_to_expression!(quad.knowns, new_coef, new_kvar)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant, new_vae::_VAE)
    JuMP.add_to_expression!(quad.variables, new_coef, new_vae)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_vae::_VAE, new_coef::_Constant)
    JuMP.add_to_expression!(quad.variables, new_coef, new_vae)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant, new_dae::_DAE)
    JuMP.add_to_expression!(quad.decisions, new_coef, new_dae)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_dae::_DAE, new_coef::_Constant)
    JuMP.add_to_expression!(quad.decisions, new_coef, new_dae)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant, new_kae::_KAE)
    JuMP.add_to_expression!(quad.knowns, new_coef, new_kae)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_kae::_KAE, new_coef::_Constant)
    JuMP.add_to_expression!(quad.knowns, new_coef, new_kae)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant, new_vqe::_VQE)
    JuMP.add_to_expression!(quad.variables, new_coef, new_vqe)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_vqe::_VQE, new_coef::_Constant)
    JuMP.add_to_expression!(quad.variables, new_coef, new_vqe)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant, new_dqe::_DQE)
    JuMP.add_to_expression!(quad.decisions, new_coef, new_dqe)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_dqe::_DQE, new_coef::_Constant)
    JuMP.add_to_expression!(quad.decisions, new_coef, new_dqe)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant, new_kqe::_KQE)
    JuMP.add_to_expression!(quad.knowns, new_coef, new_kqe)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, new_kqe::_KQE, new_coef::_Constant)
    JuMP.add_to_expression!(quad.knowns, new_coef, new_kqe)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, coef::_Constant, other::DQE)
    JuMP.add_to_expression!(quad.variables, coef, other.variables)
    JuMP.add_to_expression!(quad.decisions, coef, other.decisions)
    JuMP.add_to_expression!(quad.knowns, coef, other.knowns)
    add_cross_terms!(terms, other) = begin
        for (key, term_coef) in other
            JuMP._add_or_set!(terms, key, coef * term_coef)
        end
    end
    add_cross_terms!(quad.cross_terms, other.cross_terms)
    add_cross_terms!(quad.known_variable_terms, other.known_variable_terms)
    add_cross_terms!(quad.known_decision_terms, other.known_decision_terms)
    return quad
end

function JuMP.add_to_expression!(quad::DQE, other::DQE, coef::_Constant)
    JuMP.add_to_expression!(quad, coef, other)
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 var_1::VariableRef,
                                 var_2::Union{DecisionRef, KnownRef})
    return JuMP.add_to_expression!(quad, 1.0, var_1, var_2)
end

function JuMP.add_to_expression!(quad::DQE,
                                 var_1::VariableRef,
                                 var_2::KnownRef)
    return JuMP.add_to_expression!(quad, 1.0, var_1, var_2)
end

function JuMP.add_to_expression!(quad::DQE,
                                 var_1::DecisionRef,
                                 var_2::Union{VariableRef,DecisionRef,KnownRef})
    return JuMP.add_to_expression!(quad, 1.0, var_1, var_2)
end

function JuMP.add_to_expression!(quad::DQE,
                                 var_1::KnownRef,
                                 var_2::Union{VariableRef,DecisionRef,KnownRef})
    return JuMP.add_to_expression!(quad, 1.0, var_1, var_2)
end

function JuMP.add_to_expression!(quad::DQE,
                                 var::VariableRef,
                                 aff::_VAE)
    JuMP.add_to_expression!(quad.variables, var, aff)
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 dvar::DecisionRef,
                                 aff::_VAE)
    for (coef, term_var) in linear_terms(aff)
        key = DecisionCrossTerm(dvar, term_var)
        JuMP._add_or_set!(quad.cross_terms, key, coef)
    end
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 kvar::KnownRef,
                                 aff::_VAE)
    for (coef, term_var) in linear_terms(aff)
        key = KnownVariableCrossTerm(kvar, term_var)
        JuMP._add_or_set!(quad.known_variable_terms, key, coef)
    end
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 var::VariableRef,
                                 aff::_DAE)
    for (coef, term_dvar) in linear_terms(aff)
        key = DecisionCrossTerm(term_dvar, var)
        JuMP._add_or_set!(quad.cross_terms, key, coef)
    end
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 dvar::DecisionRef,
                                 aff::_DAE)
    JuMP.add_to_expression!(quad.decisions, dvar, aff)
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 kvar::KnownRef,
                                 aff::_DAE)
    for (coef, term_dvar) in linear_terms(aff)
        key = KnownDecisionCrossTerm(kvar, term_dvar)
        JuMP._add_or_set!(quad.known_decision_terms, key, coef)
    end
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 var::VariableRef,
                                 aff::_KAE)
    for (coef, term_kvar) in linear_terms(aff)
        key = KnownVariableCrossTerm(term_kvar, var)
        JuMP._add_or_set!(quad.known_variable_terms, key, coef)
    end
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 dvar::DecisionRef,
                                 aff::_KAE)
    for (coef, term_kvar) in linear_terms(aff)
        key = KnownDecisionCrossTerm(term_kvar, dvar)
        JuMP._add_or_set!(quad.known_decision_terms, key, coef)
    end
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 kvar::KnownRef,
                                 aff::_KAE)
    JuMP.add_to_expression!(quad.knowns, kvar, aff)
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 var::Union{VariableRef, DecisionRef, KnownRef},
                                 aff::_DAE)
    JuMP.add_to_expression!(quad, var, aff.variables)
    JuMP.add_to_expression!(quad, var, aff.decisions)
    JuMP.add_to_expression!(quad, var, aff.knowns)
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 aff::Union{_VAE, _DAE, _KAE, DAE},
                                 var::Union{VariableRef, DecisionRef, KnownRef})
    return JuMP.add_to_expression!(quad, var, aff)
end

function JuMP.add_to_expression!(quad::DQE,
                                 lhs::_VAE,
                                 rhs::_VAE)
    JuMP.add_to_expression!(quad.variables, lhs, rhs)
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 lhs::_DAE,
                                 rhs::_DAE)
    JuMP.add_to_expression!(quad.decisions, lhs, rhs)
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 lhs::_KAE,
                                 rhs::_KAE)
    JuMP.add_to_expression!(quad.knowns, lhs, rhs)
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 lhs::Union{_VAE,_DAE,_KAE},
                                 rhs::Union{_VAE,_DAE,_KAE})
    lhs_length = length(linear_terms(lhs))
    rhs_length = length(linear_terms(rhs))

    # Quadratic terms
    for (lhscoef, lhsvar) in linear_terms(lhs)
        for (rhscoef, rhsvar) in linear_terms(rhs)
            add_to_expression!(quad, lhscoef*rhscoef, lhsvar, rhsvar)
        end
    end

    # Try to preallocate space for aff
    cur = length(linear_terms(quad, typeof(rhs)))
    if !iszero(lhs.constant)
        sizehint!(quad, cur + lhs_length + rhs_length, typeof(rhs))
    end
    if !iszero(rhs.constant)
        sizehint!(quad, cur + lhs_length, typeof(lhs))
    end

    # [LHS constant] * [RHS linear terms]
    if !iszero(lhs.constant)
        c = lhs.constant
        for (rhscoef, rhsvar) in linear_terms(rhs)
            if rhs isa _VAE
                add_to_expression!(quad.variables.aff, c*rhscoef, rhsvar)
            elseif rhs isa _DAE
                add_to_expression!(quad.decisions.aff, c*rhscoef, rhsvar)
            elseif rhs isa _KAE
                add_to_expression!(quad.knowns.aff, c*rhscoef, rhsvar)
            end
        end
    end

    # [RHS constant] * [LHS linear terms]
    if !iszero(rhs.constant)
        c = rhs.constant
        for (lhscoef, lhsvar) in linear_terms(lhs)
            if lhs isa _VAE
                add_to_expression!(quad.variables.aff, c*lhscoef, lhsvar)
            elseif lhs isa _DAE
                add_to_expression!(quad.decisions.aff, c*lhscoef, lhsvar)
            elseif lhs isa _KAE
                add_to_expression!(quad.knowns.aff, c*lhscoef, lhsvar)
            end
        end
    end

    quad.variables.aff.constant += lhs.constant * rhs.constant

    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 lhs::Union{_VAE, _DAE, _KAE},
                                 rhs::DAE)
    JuMP.add_to_expression!(quad, lhs, rhs.variables)
    JuMP.add_to_expression!(quad, lhs, rhs.decisions)
    JuMP.add_to_expression!(quad, lhs, rhs.knowns)
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 lhs::DAE,
                                 rhs::Union{_VAE, _DAE, _KAE})
    JuMP.add_to_expression!(quad, rhs, lhs)
    return quad
end

function JuMP.add_to_expression!(quad::DQE,
                                 lhs::DAE,
                                 rhs::DAE)
    # Variables
    JuMP.add_to_expression!(quad, lhs.variables, rhs.variables)
    JuMP.add_to_expression!(quad, lhs.variables, rhs.decisions)
    JuMP.add_to_expression!(quad, lhs.variables, rhs.knowns)
    # Decisions
    JuMP.add_to_expression!(quad, lhs.decisions, rhs.variables)
    JuMP.add_to_expression!(quad, lhs.decisions, rhs.decisions)
    JuMP.add_to_expression!(quad, lhs.decisions, rhs.knowns)
    # Knowns
    JuMP.add_to_expression!(quad, lhs.knowns, rhs.variables)
    JuMP.add_to_expression!(quad, lhs.knowns, rhs.decisions)
    JuMP.add_to_expression!(quad, lhs.knowns, rhs.knowns)
    return quad
end

# With three factors.
function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant,
                                 new_var1::VariableRef,
                                 new_var2::VariableRef)
    JuMP.add_to_expression!(quad.variables, new_coef, new_var1, new_var2)
    return quad
end
function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant,
                                 new_var1::VariableRef,
                                 new_var2::DecisionRef)
    key = DecisionCrossTerm(new_var2, new_var1)
    JuMP._add_or_set!(quad.cross_terms, key, new_coef)
    return quad
end
function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant,
                                 new_var1::VariableRef,
                                 new_var2::KnownRef)
    key = KnownVariableCrossTerm(new_var2, new_var1)
    JuMP._add_or_set!(quad.known_variable_terms, key, new_coef)
    return quad
end
function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant,
                                 new_var1::DecisionRef,
                                 new_var2::VariableRef)
    key = DecisionCrossTerm(new_var1, new_var2)
    JuMP._add_or_set!(quad.cross_terms, key, new_coef)
    return quad
end
function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant,
                                 new_var1::DecisionRef,
                                 new_var2::DecisionRef)
    JuMP.add_to_expression!(quad.decisions, new_coef, new_var1, new_var2)
    return quad
end
function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant,
                                 new_var1::DecisionRef,
                                 new_var2::KnownRef)
    key = KnownDecisionCrossTerm(new_var2, new_var1)
    JuMP._add_or_set!(quad.known_decision_terms, key, new_coef)
    return quad
end
function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant,
                                 new_var1::KnownRef,
                                 new_var2::VariableRef)
    key = KnownVariableCrossTerm(new_var1, new_var2)
    JuMP._add_or_set!(quad.known_variable_terms, key, new_coef)
    return quad
end
function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant,
                                 new_var1::KnownRef,
                                 new_var2::DecisionRef)
    key = KnownDecisionCrossTerm(new_var1, new_var2)
    JuMP._add_or_set!(quad.known_decision_terms, key, new_coef)
    return quad
end
function JuMP.add_to_expression!(quad::DQE, new_coef::_Constant,
                                 new_var1::KnownRef,
                                 new_var2::KnownRef)
    JuMP.add_to_expression!(quad.knowns, new_coef, new_var1, new_var2)
    return quad
end


text:
"""
controls changing plane for example from transverse to saggital ...
"""
module ChangePlane
using ModernGL, ..DisplayWords, ..StructsManag, Setfield, ..PrepareWindow,   ..DataStructs , Rocket, GLFW,Dictionaries,  ..ForDisplayStructs, ..TextureManag,  ..OpenGLDisplayUtils,  ..Uniforms, Match, Parameters,DataTypesBasic

"""
In case we want to change the dimansion of scrolling so for example from transverse to coronal ...
    toBeSavedForBack - just marks weather we wat to save the info how to undo latest action
    - false if we invoke it from undoing
"""

function processKeysInfo(toScrollDatPrim::Identity{DataToScrollDims}
                    ,actor::SyncActor{Any, ActorWithOpenGlObjects}
                    ,keyInfo::KeyboardStruct
                    ,toBeSavedForBack::Bool = true ) where T
    toScrollDat= toScrollDatPrim.value

    old = actor.actor.onScrollData.dimensionToScroll

    newCalcDim= getHeightToWidthRatio(actor.actor.calcDimsStruct,toScrollDat )|>
                    getMainVerticies
     actor.actor.calcDimsStruct = newCalcDim
#In order to make the  background black  before we will render quad of possibly diffrent dimensions we will set all to invisible - and obtain black background
textSpecs = actor.actor.mainForDisplayObjects.listOfTextSpecifications

for textSpec in textSpecs
    setTextureVisibility(false,textSpec.uniforms )
end#for
basicRender(actor.actor.mainForDisplayObjects.window)

    #we need to change textures only if dimensions do not match
  #  if(actor.actor.calcDimsStruct.imageTextureWidth!=newCalcDim.imageTextureWidth  || actor.actor.calcDimsStruct.imageTextureHeight!=newCalcDim.imageTextureHeight )
        # first we need to update information about dimensions etc

        #next we need to delete all textures and create new ones

        arr = map(it->it.ID[],textSpecs)
glFinish()# make open gl ready for work

glDeleteTextures(length(arr), arr)# deleting

        #getting new
        initializeTextures(textSpecs,newCalcDim)

   # end#if

actor.actor.onScrollData.dimensionToScroll = toScrollDat.dimensionToScroll
actor.actor.onScrollData.dataToScrollDims = toScrollDat

actor.actor.onScrollData.slicesNumber = getSlicesNumber(actor.actor.onScrollData)
#getting  the slice of intrest based on last recorded mouse position

current=actor.actor.lastRecordedMousePosition[toScrollDat.dimensionToScroll]

#displaying all

singleSlDat= actor.actor.onScrollData.dataToScroll|>
(scrDat)-> map(threeDimDat->threeToTwoDimm(threeDimDat.type,Int64(current),toScrollDat.dimensionToScroll,threeDimDat ),scrDat) |>
(twoDimList)-> SingleSliceDat(listOfDataAndImageNames=twoDimList
                            ,sliceNumber=current
                            ,textToDisp = getTextForCurrentSlice(actor.actor.onScrollData, Int32(current))  )

                            # glFinish()
                            # glFlush()
                            # glClearColor(0.0, 0.0, 0.0 , 1.0)
                            # GLFW.SwapBuffers(actor.actor.mainForDisplayObjects.window)

dispObj = actor.actor.mainForDisplayObjects
#for displaying new quad - to accomodate new proportions
reactivateMainObj(dispObj.shader_program, dispObj.vbo,newCalcDim  )

glClear(GL_COLOR_BUFFER_BIT)
actor.actor.currentlyDispDat=singleSlDat = singleSlDat
updateImagesDisplayed(singleSlDat
                    ,actor.actor.mainForDisplayObjects
                    ,actor.actor.textDispObj
                    ,newCalcDim
                    ,actor.actor.valueForMasToSet      )

     #saving information about current slice for future reference
actor.actor.currentDisplayedSlice = current
# to enbling getting back
if(toBeSavedForBack)
    addToforUndoVector(actor, ()-> processKeysInfo( Option(old),actor, keyInfo,false ))
end#if

end#processKeysInfo
end#ChangePlane


text:
!(dirname(@__FILE__) in LOAD_PATH) && push!(LOAD_PATH, dirname(@__FILE__))
using DSP, Test, Random, FilterTestHelpers
#
# filt with different filter forms
#
@testset "filt" begin
    x = randn(100)

    for n = 1:6
        for proto in (Butterworth(n), Chebyshev1(n, 1), Chebyshev2(n, 1), Elliptic(n, 0.5, 2))
            zpk = digitalfilter(Lowpass(0.2), proto)
            tf = convert(PolynomialRatio, zpk)
            if n <= 2
                bq = convert(Biquad, zpk)
            else
                @test_throws ArgumentError convert(Biquad, zpk)
            end
            sos = convert(SecondOrderSections, zpk)

            res = filt(sos, x)

            # Test with filt with tf/sos
            tfres = filt(tf, x)
            @test res ≈ tfres
            @test res ≈ filt!(similar(x), sos, x)
            @test res ≈ filt!(similar(x), tf, x)

            # For <= 2 poles, test with biquads
            if n <= 2
                @test res ≈ filt(bq, x)
                @test res ≈ filt!(similar(x), bq, x)
                f = DF2TFilter(bq)
                @test tfres ≈ [filt(f, x[1:50]); filt(f, x[51:end])]
            end

            # Test that filt with zpk converts
            @test res ≈ filt(zpk, x)
            @test res ≈ filt!(similar(x), zpk, x)

            # Test with DF2TFilter
            f = DF2TFilter(sos)
            @test res ≈ [filt(f, x[1:50]); filt(f, x[51:end])]
            f = DF2TFilter(tf)
            @test tfres ≈ [filt(f, x[1:50]); filt(f, x[51:end])]
            f = DF2TFilter(zpk)
            @test res ≈ [filt(f, x[1:50]); filt(f, x[51:end])]
        end
    end

    # Test simple scaling with DF2TFilter
    @test filt(DF2TFilter(PolynomialRatio([3.7], [4.2])), x) == x * (3.7/4.2)
end

#
# filtfilt
#

@testset "filt_stepstate" begin
    ##############
    #
    # Filter initial conditions
    # Python example 1 - http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter_zi.html
    #
    # b, a = butter(5, 0.25)
    # zi = lfilter_zi(b, a)
    #
    ##############
    zi_python = [ 0.99672078, -1.49409147,  1.28412268, -0.45244173,  0.07559489]

    b = [ 0.00327922,  0.01639608,  0.03279216,  0.03279216,  0.01639608,  0.00327922]
    a = [ 1.        , -2.47441617,  2.81100631, -1.70377224,  0.54443269, -0.07231567]

    @test ≈(zi_python, DSP.Filters.filt_stepstate(b, a), atol=1e-7)

    ##############
    #
    # Filter initial conditions
    # Matlab - Random filter coefficients
    #
    # filtfilt([0.222, 0.43, 0.712], [1, 0.33, 0.22], x)
    # Then used breakpoints to extract zi
    #
    ##############

    zi_matlab = [0.6580, 0.5184]

    b = [0.222, 0.43, 0.712]
    a = [1, 0.33, 0.22]

    @test zi_matlab ≈ DSP.Filters.filt_stepstate(b, a)

    ##############
    #
    # Filter initial conditions
    # Python - non 1 first coeffecient
    #
    # b = array([ 0.00327922,  0.01639608,  0.03279216,  0.03279216,  0.01639608,  0.00327922])
    # a = array([ 1.1       , -2.47441617,  2.81100631, -1.70377224,  0.54443269, -0.07231567])
    # zi = lfilter_zi(b, a)
    #
    ##############

    zi_python = [0.55996501, -0.72343165,  0.68312446, -0.2220676 ,  0.04030775]

    b = [ 0.00327922,  0.01639608,  0.03279216,  0.03279216,  0.01639608,  0.00327922]
    a = [ 1.1       , -2.47441617,  2.81100631, -1.70377224,  0.54443269, -0.07231567]

    @test ≈(zi_python, DSP.Filters.filt_stepstate(b, a), atol=1e-7)
end

##############
#
# Filter in place check (with initial conditions)
#
# x = '/Users/rluke/.julia/v0.3/DSP/test/data/spectrogram_x.txt';
# filtcheck  = filter([0.4, 1], [0.9, 0.6], x, 0.4750)
# dlmwrite('filt_check.txt',[filtcheck], 'delimiter', '\t', 'precision', '%.12f')
#
##############
@testset "filt! with init. cond." begin
    matlab_filt  = readdlm(joinpath(dirname(@__FILE__), "data", "filt_check.txt"),'\t')

    a = [0.9, 0.6]
    b = [0.4, 1]
    z = [0.4750]
    x  = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t')
    DSP.Filters.filt!(vec(x), b, a, vec(x), z)

    @test matlab_filt ≈ x
end

#######################################
#
# Test 1d filtfilt against matlab results
#
# x = '/Users/rluke/.julia/v0.3/DSP/test/data/spectrogram_x.txt'; x = textread(x);
# b = [ 0.00327922,  0.01639608,  0.03279216,  0.03279216,  0.01639608,  0.00327922];
# a = [ 1.        , -2.47441617,  2.81100631, -1.70377224,  0.54443269, -0.07231567];
# x2 = filtfilt(b, a, x);
# dlmwrite('filtfilt_output.txt',[x2], 'delimiter', '\t', 'precision', '%.12f')
#
#######################################
@testset "filtfilt 1D" begin
    x2_matlab = readdlm(joinpath(dirname(@__FILE__), "data", "filtfilt_output.txt"),'\t')

    b = [ 0.00327922,  0.01639608,  0.03279216,  0.03279216,  0.01639608,  0.00327922]
    a = [ 1.        , -2.47441617,  2.81100631, -1.70377224,  0.54443269, -0.07231567]
    x  = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t')

    @test x2_matlab ≈ filtfilt(b, a, x)
end

# Make sure above doesn't crash for real coeffs & complex data.
@testset "filtfilt 1D Complex" begin
    x2_matlab = readdlm(joinpath(dirname(@__FILE__), "data", "filtfilt_output.txt"),'\t')

    b = [ 0.00327922,  0.01639608,  0.03279216,  0.03279216,  0.01639608,  0.00327922]
    a = [ 1.        , -2.47441617,  2.81100631, -1.70377224,  0.54443269, -0.07231567]
    x  = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t')

    y = x .+ 1im .* randn(size(x, 1))

    @test x2_matlab ≈ real.(filtfilt(b, a, y))
end

#######################################
#
# Test 2d filtfilt against matlab results
#

#=x = '/Users/rluke/.julia/v0.3/DSP/test/data/spectrogram_x.txt'; x = textread(x);=#

#=x = repmat(x, 1, 3);=#
#=x(:,2) = circshift(x(:,2), 64);=#
#=x(:,3) = circshift(x(:,3), 128);=#

#=b = [ 0.00327922,  0.01639608,  0.03279216,  0.03279216,  0.01639608,  0.00327922];=#
#=a = [ 1.        , -2.47441617,  2.81100631, -1.70377224,  0.54443269, -0.07231567];=#
#=x2 = filtfilt(b, a, x);=#
#=dlmwrite('filtfilt_output_2d.txt',[x2], 'delimiter', '\t', 'precision', '%.12f')=#

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

@testset "filtfilt 2D" begin
    b = [ 0.00327922,  0.01639608,  0.03279216,  0.03279216,  0.01639608,  0.00327922]
    a = [ 1.        , -2.47441617,  2.81100631, -1.70377224,  0.54443269, -0.07231567]

    x2_output = readdlm(joinpath(dirname(@__FILE__), "data", "filtfilt_output_2d.txt"),'\t')

    x  = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t')
    x  = repeat(x, outer=(1, 3))
    x[:,2] = circshift(x[:,2], 64)
    x[:,3] = circshift(x[:,3], 128)

    @test x2_output ≈ filtfilt(b, a, x)

    #######################################
    #
    # Test 2d filtfilt with filter type
    #
    #
    #######################################

    f = PolynomialRatio(b, a)

    # Use 2d data from last test
    @test x2_output ≈ filtfilt(f, x)
end

#######################################
#
# Test 2d filtfilt with SecondOrderSections
#
# Our implementation differs from MATLAB, but should match the
# PolynomialRatio provided the filter order is even. (Otherwise
# the extrapolation will differ slightly.)
#######################################
@testset "filtfilt SOS" begin
    x  = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t')

    f = DSP.digitalfilter(DSP.Lowpass(0.2), DSP.Butterworth(4))
    @test filtfilt(convert(SecondOrderSections, f), x) ≈ filtfilt(convert(PolynomialRatio, f), x)

    f = DSP.digitalfilter(DSP.Highpass(0.1), DSP.Butterworth(6))
    @test filtfilt(convert(SecondOrderSections, f), x) ≈ filtfilt(convert(PolynomialRatio, f), x)

    f = DSP.digitalfilter(DSP.Bandpass(0.1, 0.3), DSP.Butterworth(2))
    @test filtfilt(convert(SecondOrderSections, f), x) ≈ filtfilt(convert(PolynomialRatio, f), x)
end

#
# fftfilt/filt
#

@testset "fftfilt $xlen/$blen" for xlen in 2 .^ (7:18) .- 1, blen in 2 .^ (1:6) .- 1
    b = randn(blen)
    for x in (rand(xlen), rand(xlen, 2))
        out = similar(x)
        filtres = filt(b, [1.0], x)
        fftres = fftfilt(b, x)
        firres = filt(b, x)
        td_res = tdfilt(b, x)

        @test filtres ≈ fftres
        @test filtres ≈ firres
        @test filtres ≈ td_res

        fftfilt!(out, b, x)
        @test filtres ≈ out # test output of fftfilt!

        tdfilt!(out, b, x)
        @test filtres ≈ out # test output of tdfilt!
    end
end

# fir_filtfilt
@testset "fir_filtfilt" begin
    b = randn(10)
    for x in (randn(100), randn(100, 2))
        @test DSP.Filters.filtfilt(b, x) ≈ DSP.Filters.iir_filtfilt(b, [1.0], x)
        @test DSP.Filters.filtfilt(b, [2.0], x) ≈ DSP.Filters.iir_filtfilt(b, [2.0], x)
    end
end


text:
include("BEMTypes.jl")

function assemble_square_model()
    # Assemble mesh
    dim = Dimension2()

    vertices = Point[]
    push!(vertices, Point(0.0, 0.0))
    push!(vertices, Point(1.0, 0.0))
    push!(vertices, Point(1.0, 1.0))
    push!(vertices, Point(0.0, 1.0))

    # Initialize empty nodes array. Nodes are added during discretization.
    nodes = Point[]
    # Initialize empty normal array. Normal vectors are added during
    # discretization.
    normals = Point[]
    geometry = MeshGeometry(vertices, nodes, normals)

    # Mesh connectivity dictionary
    connectivity = Dict{AbstractString, MeshConnectivity}()
    # Edge to vertex connectivity
    indices = [1, 2, 2, 3, 3, 4, 4, 1]
    offsets = [1, 3, 5, 7, 9]
    connectivity["E2V"] = MeshConnectivity(indices, offsets)
    topology = MeshTopology(dim, connectivity)

    mesh = Mesh(topology, geometry)

    # Assemble boundary condition
    ice(x,y) = 0.0
    fire(x,y) = 100.0
    insulated(x,y) = 0.0
    ice_boundary = Dirichlet(ice)
    fire_boundary = Dirichlet(fire)
    insulated_boundary = Neumann(insulated)
    bcs = BoundaryCondition[ice_boundary, fire_boundary, insulated_boundary]
    markers = [1, 3, 2, 3]
    mesh_bc = MeshBoundaryCondition(bcs, markers)

    return mesh, mesh_bc
end

function assemble_cube_model()
    # Assemble mesh
    dim = Dimension3()

    vertices = Point[]
    push!(vertices, Point(0.0, 0.0, 0.0))
    push!(vertices, Point(1.0, 0.0, 0.0))
    push!(vertices, Point(1.0, 0.0, 1.0))
    push!(vertices, Point(0.0, 0.0, 1.0))
    push!(vertices, Point(1.0, 1.0, 0.0))
    push!(vertices, Point(1.0, 1.0, 1.0))
    push!(vertices, Point(0.0, 1.0, 0.0))
    push!(vertices, Point(0.0, 1.0, 0.0))

    # Initialize empty nodes array. Nodes are added during discretization.
    nodes = Point[]
    # Initialize empty normal array. Normal vectors are added during
    # discretization.
    normals = Point[]
    geometry = MeshGeometry(vertices, nodes, normals)

    # Mesh connectivity dictionary
    connectivity = Dict{AbstractString, MeshConnectivity}()
    # Face to vertex connectivity
    indices = [1, 3, 4, 1, 2, 3, 2, 5, 3, 5, 6, 3, 5, 7, 6, 8, 7, 5, 8, 4, 7, 1,
               8, 4, 4, 3, 7, 3, 6, 7, 1, 2, 5, 1, 5, 8]
    offsets = [1, 4, 7, 10, 13, 16, 19, 22, 25, 28, 31, 34, 37]
    F2V = MeshConnectivity(indices, offsets)
    connectivity["F2V"] = F2V
    topology = MeshTopology(dim, connectivity)

    mesh = Mesh(topology, geometry)

    # Assemble boundary condition
    ice(x,y,z) = 0.0
    fire(x,y,z) = 100.0
    insulated(x,y,z) = 0.0
    ice_boundary = Dirichlet(ice)
    fire_boundary = Dirichlet(fire)
    insulated_boundary = Neumann(insulated)
    bcs = BoundaryCondition[ice_boundary, fire_boundary, insulated_boundary]
    markers = [3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 1, 1]
    mesh_bc = MeshBoundaryCondition(bcs, markers)

    return mesh, mesh_bc
end


text:
## Test grenander and related methods ##
module Test_grenander

using MultipleTesting
using Base.Test

## Grenander ##
println(" ** ", "grenander")
pv = [0.2, 0.4, 0.6, 0.8, 1.0]
p, f, F = MultipleTesting.grenander(pv)
@test_approx_eq p pv
@test_approx_eq f ones(p)
@test_approx_eq F pv

pv = [0.1, 0.3, 0.5, 0.7, 0.9]
p, f, F = MultipleTesting.grenander(pv)
@test_approx_eq p pv
@test_approx_eq f ones(p)
@test_approx_eq F pv+0.1

pv = [0.1, 0.1, 0.3, 0.3, 0.5, 0.5, 0.7, 0.7, 0.9, 0.9]
p, f, F = MultipleTesting.grenander(pv)
@test_approx_eq p unique(pv)
@test_approx_eq f ones(p)
@test_approx_eq F p+0.1

pv = [0.1, 0.2, 0.5, 0.6, 0.9]
p, f, F = MultipleTesting.grenander(pv)
@test_approx_eq p pv
@test_approx_eq f [2.0, 1.0, 1.0, 2/3, 2/3]
@test_approx_eq F [0.2, 0.4, 0.7, 0.8, 1.0]

## isotonic regression ##
println(" ** ", "isotonic regression")

#pooled adjacent violators example page 10 robertson
isotonicregression = MultipleTesting.isotonicregression
@test_approx_eq isotonicregression([22.5; 23.333; 20.833; 24.25], [3.0;3.0;3.0;2.0]) [22.222; 22.222; 22.222; 24.25]
#if input already ordered, then output should be the same
@test_approx_eq isotonicregression([1.0; 2.0; 3.0]) [1.0; 2.0; 3.0]
@test_approx_eq isotonicregression([1., 41., 51., 1., 2., 5., 24.], [1., 2., 3., 4., 5., 6., 7.]) [1.0, 13.95, 13.95, 13.95, 13.95, 13.95, 24]

#pooled adjacent violators example page 10 robertson
isotonic_regression = MultipleTesting.isotonic_regression
@test_approx_eq isotonic_regression([22.5; 23.333; 20.833; 24.25], [3.0;3.0;3.0;2.0]) [22.222; 22.222; 22.222; 24.25]
#if input already ordered, then output should be the same
@test_approx_eq isotonic_regression([1.0; 2.0; 3.0]) [1.0; 2.0; 3.0]
@test_approx_eq isotonic_regression([1., 41., 51., 1., 2., 5., 24.], [1., 2., 3., 4., 5., 6., 7.]) [1.0, 13.95, 13.95, 13.95, 13.95, 13.95, 24]

# single value or empty vector remains unchanged
r = rand(1)
@test MultipleTesting.isotonic_regression(r) == r
r = Vector{Float64}()
@test MultipleTesting.isotonic_regression(r) == r

r = rand(10)
@test MultipleTesting.isotonic_regression(r, ones(r)) == MultipleTesting.isotonic_regression(r)
@test_throws DimensionMismatch MultipleTesting.isotonic_regression(rand(10), ones(5))

end


text:
push!(LOAD_PATH, ".")

module Individual

include("utility.jl")

abstract type AbstractIndividual{GeneType} <: AbstractArray{GeneType, 1} end

@define Individual begin
  chromosome::Vector{GeneType}
  num_genes::Integer
  bounds::Vector{Tuple{GeneType, GeneType}}
end

function build(ng::Integer, b::Vector{Tuple{GeneType, GeneType}}) where {GeneType}
  return Vector{GeneType}(undef, ng), ng, b
end

function build(ind::AbstractIndividual)
  return deepcopy(ind.chromosome), deepcopy(ind.num_genes), deepcopy(ind.bounds)
end

function generateRandom!() end

function getGeneValue(this::AbstractIndividual, pos::Integer)
  return this.chromosome[pos]
end

function getNumGenes(this::AbstractIndividual)
  return this.num_genes
end

function getGeneBounds(this::AbstractIndividual, pos::Integer)
    return this.bounds[pos]
end

function setGeneValue!(this::AbstractIndividual{GeneType}, pos::Integer, value::GeneType) where {GeneType}
    return this.chromosome[pos] = value
end

function setGeneBounds!(this::AbstractIndividual{GeneType}, pos::Integer, bound::Vector{Tuple{GeneType, GeneType}}) where {GeneType}
    return this.bounds[pos] = bound
end

function Base.show(io::IO, this::AbstractIndividual)
  println(io, "$(this.num_genes) genes.")
  println(io, "Chromosome: ", this.chromosome)
  println(io, "Bounds: ", this.bounds)
end

function Base.size(this::AbstractIndividual)
  return (getNumGenes(this), )
end

function Base.similar(this::AbstractIndividual, ::Type{GeneType}, sz::Integer) where {GeneType}
  return Vector{GeneType}(undef, sz)
end

function Base.IndexStyle(::Type{<:AbstractIndividual})
  return IndexLinear()
end

function Base.getindex(this::AbstractIndividual, pos::Integer)
  return this.chromosome[pos]
end

function Base.setindex!(this::AbstractIndividual{GeneType}, value::GeneType, pos::Integer) where {GeneType}
  return this.chromosome[pos] = value
end

function Base.firstindex(this::AbstractIndividual)
  return 1
end

function Base.lastindex(this::AbstractIndividual)
  return getNumGenes(this)
end

end

max_stars_repo_name: Markhyz/Julia-Genetic-Algorithm

max_issues_repo_name: Markhyz/Julia-Genetic-Algorithm

max_forks_repo_name: Markhyz/Julia-Genetic-Algorithm


text:
# tolerance is ||alpha^k - alpha^k+1||_infty
mutable struct FIBSolver <: Solver
    max_iterations::Int64
    tolerance::Float64
    verbose::Bool
end
function FIBSolver(;max_iterations::Int64=100, tolerance::Float64=1e-3, verbose::Bool=false)
    return FIBSolver(max_iterations, tolerance, verbose)
end

function solve(solver::FIBSolver, pomdp::POMDP; kwargs...)
    if !isempty(kwargs)
        @warn("Keyword args for solve(::FIBSolver, ::MDP) are no longer supported. For verbose output, use the verbose option in the FIBSolver")
    end

    state_list = ordered_states(pomdp)
    obs_list = observations(pomdp)
    action_list = ordered_actions(pomdp)

    ns = length(state_list)
    na = length(action_list)

    alphas = zeros(ns,na)
    old_alphas = zeros(ns,na)

    for i = 1:solver.max_iterations

        # copyto!(dest, src)
        copyto!(old_alphas, alphas)

        residual = 0.0

        for (ai, a) in enumerate(action_list)
            for (si, s) in enumerate(state_list)

                sp_dist = transition(pomdp, s, a)

                r = 0.0
                for (sp, p_sp) in weighted_iterator(sp_dist)
                    r += p_sp*reward(pomdp, s, a, sp)
                end

                # Sum_o max_a' Sum_s' O(o | s',a) T(s'|s,a) alpha_a^k(s')
                o_sum = 0.0
                for o in obs_list

                    # take maximum over ap
                    ap_sum = -Inf
                    for (api, ap) in enumerate(action_list)

                        # Sum_s' O(o | s',a) T(s'|s,a) alpha_a^k(s')
                        temp_ap_sum = 0.0
                        for (sp, p_sp) in weighted_iterator(sp_dist)
                            o_dist = observation(pomdp, s, a, sp)
                            p_o = pdf(o_dist, o)
                            spi = stateindex(pomdp, sp)

                            temp_ap_sum += p_o * p_sp * old_alphas[spi,api]
                        end
                        ap_sum = max(temp_ap_sum, ap_sum)
                    end

                    o_sum += ap_sum
                end

                alphas[si, ai] = r + discount(pomdp) * o_sum

                alpha_diff = abs(alphas[si, ai] - old_alphas[si, ai])
                residual = max(alpha_diff, residual)
            end
        end

        solver.verbose ? @printf("[Iteration %-4d] residual: %10.3G \n", i, residual) : nothing
        residual < solver.tolerance ? break : nothing
    end

    return AlphaVectorPolicy(pomdp, alphas, action_list)
end


text:
## Transforms for MalmquistTakenaka Space
hasfasttransform(::MalmquistTakenaka) = true
# note that one of these points is always infinity
points(S::MalmquistTakenaka{T},n::Int) where T = real(S.λ) .+ imag(S.λ)*tan.(range(zero(real(T)),length=n,step=π/n))

function plan_transform(S::MalmquistTakenaka,vals::AbstractVector)
    weights = sqrt(π/abs(imag(S.λ))).*(conj(S.λ) .- points(S,length(vals)))
    TransformPlan(S,(weights,FastTransforms.plan_fft(complex(vals))),Val{false})
end
function plan_itransform(S::MalmquistTakenaka,cfs::AbstractVector)
    weights = sqrt(abs(imag(S.λ))/π) ./ (conj(S.λ) .- points(S,length(cfs)))
    ITransformPlan(S,(weights,FastTransforms.plan_bfft(complex(cfs))),Val{false})
end
function *(P::TransformPlan{T,S,false},vals::AbstractVector) where {T,S<:MalmquistTakenaka}
    weights,fftplan = P.plan
    _mt_reorder_scale_and_phase_shift(fftplan*(weights.*complex(vals)))
end
function *(P::ITransformPlan{T,S,false},cfs::AbstractVector) where {T,S<:MalmquistTakenaka}
    weights,ifftplan = P.plan
    vals = weights.*(ifftplan*(_imt_reorder_and_phase_shift(complex(cfs))))
    if (iseven(length(vals))) vals[div(length(vals),2)+1] = zero(eltype(vals)) end
    return vals
end

## Transforms for Weideman Space
hasfasttransform(::Weideman) = true
# note that one of these points is always infinity
points(S::Weideman{T},n::Int) where T = real(S.λ) .+ imag(S.λ)*tan.(range(zero(real(T)) ,length=n,step=π/n))

function plan_transform(S::Weideman,vals::AbstractVector)
    TransformPlan(S,FastTransforms.plan_fft(complex(vals)),Val{false})
end
function plan_itransform(S::Weideman,cfs::AbstractVector)
    ITransformPlan(S,FastTransforms.plan_bfft(complex(cfs)),Val{false})
end
function *(P::TransformPlan{T,S,false},vals::AbstractVector) where {T,S<:Weideman}
    fftplan = P.plan
    _mt_reorder_scale_and_phase_shift(fftplan*((-one(T)im)*complex(vals)))
end
function *(P::ITransformPlan{T,S,false},cfs::AbstractVector) where {T,S<:Weideman}
    ifftplan = P.plan
    (one(T)im)*ifftplan*(_imt_reorder_and_phase_shift(complex(cfs)))
end

# Deals with the fact that the positively and negatively indexed
# coefficients need to be interlaced
function _mt_reorder_scale_and_phase_shift(cfs::Vector{T}) where T
    n = length(cfs)
    tmp = FastTransforms.fftshift(cfs)
    tmp .*= (-one(T)*im).^(range(-floor(n/2)-one(T),length=n))/n
    newcfs = copy(cfs)
    newcfs[1:2:n] = tmp[div(n,2)+1:n]
    newcfs[2:2:n] = tmp[div(n,2):-1:1]
    newcfs
end

# Deals with the fact that the positively and negatively indexed
# coefficients need to be interlaced
function _imt_reorder_and_phase_shift(cfs::Vector{T}) where T
    n = length(cfs)
    newcfs = copy(cfs)
    newcfs[div(n,2)+1:n] = cfs[1:2:n]
    newcfs[div(n,2):-1:1] = cfs[2:2:n]
    newcfs .*= (one(T)*im).^(range(-floor(n/2)-one(T),length=n))
    FastTransforms.ifftshift(newcfs)
end


text:
# Testing plot examples
using Plots
using SmoothLivePlot
using Test
include("testPlotFunctions.jl")
gr(show = true)

# Test macro
macro no_error(ex)
    quote
        try
            $(esc(ex))
            true
        catch
            false
        end
    end
end

function main()
    @testset "All plots" begin
        @testset "Modify Y" begin
            @test @no_error testModifyY()
        end
        @testset "Modify X" begin
            @test @no_error testModifyX()
        end
        @testset "Modify X+Y" begin
            @test @no_error testModifyXY()
        end
        @testset "Modify Z" begin
            @test @no_error testModifyZ()
        end
        @testset "Modify X+Text" begin
            @test @no_error testModifyXText()
        end
        @testset "Modify X+Colour" begin
            @test @no_error testModifyXColour()
        end
        @testset "Add to X+Y" begin
            @test @no_error testAddXY()
        end
        @testset "Add to X+Y+Z" begin
            @test @no_error testAddXYZ()
        end
    end
end

main()


text:
# from https://github.com/chengchingwen/Transformers.jl/tree/master/src/fix

for typ ∈ atomictypes
    lt = llvmtypes[typ]
    ilt = llvmtypes[inttype(typ)]
    # Note: atomic_cas! succeeded (i.e. it stored "new") if and only if the result is "cmp"
    if typ <: Integer
        @eval atomic_cas!(x::Ptr{$typ}, cmp::$typ, new::$typ) =
            llvmcall($"""
                     %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                     %rs = cmpxchg $lt* %ptr, $lt %1, $lt %2 acq_rel acquire
                     %rv = extractvalue { $lt, i1 } %rs, 0
                     ret $lt %rv
                     """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                     x, cmp, new)
    else
        @eval atomic_cas!(x::Ptr{$typ}, cmp::$typ, new::$typ) =
            llvmcall($"""
                     %iptr = inttoptr i$WORD_SIZE %0 to $ilt*
                     %icmp = bitcast $lt %1 to $ilt
                     %inew = bitcast $lt %2 to $ilt
                     %irs = cmpxchg $ilt* %iptr, $ilt %icmp, $ilt %inew acq_rel acquire
                     %irv = extractvalue { $ilt, i1 } %irs, 0
                     %rv = bitcast $ilt %irv to $lt
                     ret $lt %rv
                     """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                     x, cmp, new)
    end

    arithmetic_ops = [:add, :sub]
    for rmwop in [arithmetic_ops..., :xchg, :and, :nand, :or, :xor, :max, :min]
        rmw = string(rmwop)
        fn = Symbol("atomic_", rmw, "!")
        if (rmw == "max" || rmw == "min") && typ <: Unsigned
            # LLVM distinguishes signedness in the operation, not the integer type.
            rmw = "u" * rmw
        end
        if rmwop in arithmetic_ops && !(typ <: ArithmeticTypes) continue end
        if typ <: Integer
            @eval $fn(x::Ptr{$typ}, v::$typ) =
                llvmcall($"""
                         %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                         %rv = atomicrmw $rmw $lt* %ptr, $lt %1 acq_rel
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, x, v)
        else
            rmwop == :xchg || continue
            @eval $fn(x::Ptr{$typ}, v::$typ) =
                llvmcall($"""
                         %iptr = inttoptr i$WORD_SIZE %0 to $ilt*
                         %ival = bitcast $lt %1 to $ilt
                         %irv = atomicrmw $rmw $ilt* %iptr, $ilt %ival acq_rel
                         %rv = bitcast $ilt %irv to $lt
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, x, v)
        end
    end
end

const opnames = Dict{Symbol, Symbol}(:+ => :add, :- => :sub, :* => :mul, :/ => :div)
for op in [:+, :-, :max, :min, :*, :/]
    opname = get(opnames, op, op)
    for typ ∈ atomictypes
        typ == Bool && continue
        !(op == :/ || op == :*) && typ <: Integer && continue
        @eval function $(Symbol("atomic_", opname, "!"))(var::Ptr{T}, val::T) where T<: $typ
            IT = inttype(T)
            old = unsafe_load(var)
            while true
                new = $op(old, val)
                cmp = old
                old = atomic_cas!(var, cmp, new)
                reinterpret(IT, old) == reinterpret(IT, cmp) && return new
                # Temporary solution before we have gc transition support in codegen.
                ccall(:jl_gc_safepoint, Cvoid, ())
            end
        end
    end
end

for op = [:add, :sub, :max, :min, :mul, :div]
    @eval function $(Symbol("scatter_", op, "!"))(ys::Matrix{T}, us::Array{T}, xs::Array{Int}) where T
        l = length(xs)
        s = size(ys, 1)

        Threads.@threads for num = 1:l*s
            li = (num -1) ÷ s + 1
            i = (num - 1) % s + 1
            @inbounds ind = Tuple(CartesianIndices(xs)[li])
            @inbounds $(Symbol("atomic_", op, "!"))(
                pointer(ys,
                        Base._to_linear_index(ys,
                                              i,
                                              xs[li]
                                              )
                        ),
                us[i, ind...]
            )
        end

        return ys
    end

    @eval function $(Symbol("scatter_", op, "!"))(ys::Array{T}, us::Array{T}, xs::Array{<:Tuple}) where T
        l = length(xs)
        s = size(ys, 1)

        Threads.@threads for num = 1:l*s
            li = (num -1) ÷ s + 1
            i = (num - 1) % s + 1
            @inbounds ind = Tuple(CartesianIndices(xs)[li])
            @inbounds $(Symbol("atomic_", op, "!"))(
                pointer(ys,
                        Base._to_linear_index(ys,
                                              i,
                                              xs[li]...
                                              )
                        ),
                us[i, ind...]
            )
        end

        return ys
    end
end

function scatter_mean!(ys::Matrix{T}, us::Array{T}, xs::Array{Int}) where T
    l = length(xs)
    s = size(ys, 1)
    N = size(ys, 2)
    Ns = [sum(xs .== i) for i = 1:N]
    ys_ = fill!(similar(ys), 0)

    Threads.@threads for num = 1:l*s
        li = (num -1) ÷ s + 1
        i = (num - 1) % s + 1
        @inbounds ind = Tuple(CartesianIndices(xs)[li])
        @inbounds atomic_add!(
            pointer(ys_, Base._to_linear_index(ys_, i, xs[li])),
            us[i, ind...]
        )
    end

    Threads.@threads for i = 1:N
        for j = 1:s
            @inbounds atomic_div!(
                pointer(ys_, Base._to_linear_index(ys_, j, i)),
                T(Ns[i])
            )
            @inbounds atomic_add!(
                pointer(ys, Base._to_linear_index(ys, j, i)),
                ys_[j,i]
            )
        end
    end

    return ys
end

# function scatter_mean!(ys::Array{T}, us::Array{T}, xs::Array{<:Tuple}) where T
#     l = length(xs)
#     s = size(ys, 1)
#     N = length(Set(xs))
#     Ns = [sum(xs .== i) for i = 1:N]
#     ys_ = fill!(similar(ys), 0)
#
#     Threads.@threads for num = 1:l*s
#         li = (num -1) ÷ s + 1
#         i = (num - 1) % s + 1
#         @inbounds ind = Tuple(CartesianIndices(xs)[li])
#         @inbounds atomic_add!(
#             pointer(ys_, Base._to_linear_index(ys_, i, xs[li]...)),
#             us[i, ind...]
#         )
#     end
#
#     Threads.@threads for i = 1:N
#         for j = 1:s
#             @inbounds atomic_div!(
#                 pointer(ys_, Base._to_linear_index(ys_, j, i)),
#                 T(Ns[i])
#             )
#             @inbounds atomic_add!(
#                 pointer(ys, Base._to_linear_index(ys, j, i)),
#                 ys_[j,i]
#             )
#         end
#     end
#
#     return ys
# end

scatter!(op::Symbol, ys::Array, us::Array, xs::Array) = scatter!(Val(op), ys, us, xs)
scatter!(::Val{:add}, ys::Array, us::Array, xs::Array) = scatter_add!(ys, us, xs)
scatter!(::Val{:sub}, ys::Array, us::Array, xs::Array) = scatter_sub!(ys, us, xs)
scatter!(::Val{:mul}, ys::Array, us::Array, xs::Array) = scatter_mul!(ys, us, xs)
scatter!(::Val{:div}, ys::Array, us::Array, xs::Array) = scatter_div!(ys, us, xs)
scatter!(::Val{:max}, ys::Array, us::Array, xs::Array) = scatter_max!(ys, us, xs)
scatter!(::Val{:min}, ys::Array, us::Array, xs::Array) = scatter_min!(ys, us, xs)
scatter!(::Val{:mean}, ys::Array, us::Array, xs::Array) = scatter_mean!(ys, us, xs)


text:
# # This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).

function cluster_detector_hits(
    detno::AbstractVector{<:Integer},
    edep::AbstractVector{<:RealQuantity},
    pos::AbstractVector{<:StaticVector{3,PT}},
    cluster_radius::RealQuantity
) where {TT<:RealQuantity,PT <: RealQuantity}
    Table = TypedTables.Table
    unsorted = Table(detno = detno, edep = edep, pos = pos)
    sorting_idxs = sortperm(unsorted.detno)
    sorted = unsorted[sorting_idxs]
    grouped = Table(consgroupedview(sorted.detno, Tables.columns(sorted)))

    r_detno = similar(detno, 0)
    r_edep = similar(edep, 0)
    r_pos = similar(pos, 0)

    posunit = unit(PT)
    ustripped_cradius = ustrip(uconvert(posunit, cluster_radius))

    for d_hits_nt in grouped
        d_hits = Table(d_hits_nt)
        d_detno = first(d_hits.detno)
        @assert all(isequal(d_detno), d_hits.detno)
        if length(d_hits) > 3
            clusters = Clustering.dbscan(ustrip.(flatview(d_hits.pos)), ustripped_cradius, leafsize = 20, min_neighbors = 1, min_cluster_size = 1)
            for c in clusters
                idxs = vcat(c.boundary_indices, c.core_indices)
                @assert length(idxs) == c.size
                c_hits = view(d_hits, idxs)
                global c_hits = c_hits

                push!(r_detno, d_detno)
                push!(r_edep, sum(c_hits.edep))
                esum = ustrip(r_edep[end])
                inv_e_sum = inv(esum)
                weights = Weights(ustrip.(c_hits.edep), esum) .* inv_e_sum
                push!(r_pos, sum(c_hits.pos .* weights))
            end
        else
            append!(r_detno, d_hits.detno)
            append!(r_edep, d_hits.edep)
            append!(r_pos, d_hits.pos)
        end
    end

    (detno = r_detno, edep = r_edep, pos = r_pos)
end

function cluster_detector_hits(table::TypedTables.Table, cluster_radius::RealQuantity)
    @assert :pos in TypedTables.columnnames(table) "Table has no column `pos`"
    @assert :edep in TypedTables.columnnames(table) "Table has no column `edep`"
    @assert :detno in TypedTables.columnnames(table) "Table has no column `detno`"
    clustered_nt = map(
        evt -> cluster_detector_hits(evt.detno, evt.edep, evt.pos, cluster_radius),
        table
    )
    TypedTables.Table(merge(
        TypedTables.columns(table),
        map(
            VectorOfVectors,
            TypedTables.columns(clustered_nt)
        )
    ))
end