text:
"""
    MIDIFile <: Any
Type representing a file of MIDI data.

## Fields
* `format::UInt16` : The format of the file. Can be 0, 1 or 2.
* `tpq::Int16` : The time division of the track, ticks-per-quarter-note.
* `tracks::Array{MIDITrack, 1}` : The array of contained tracks.
"""
mutable struct MIDIFile
    format::UInt16 # Can be 0, 1 or 2
    tpq::Int16 # The time division of the track. Ticks per quarter note
    tracks::Vector{MIDITrack}
end
# Pretty print
function Base.show(io::IO, midi::MIDIFile)
	tnames = tracknames(midi)
	s = "MIDIFile (format=$(Int(midi.format)), tpq=$(midi.tpq)) "
	if any(!isequal(NOTRACKNAME), tnames) # we have tracknames
		s *= "with tracks:\n"
		for t in tnames
			s *= " "*t*"\n"
		end
	else # some track doesn't have a name
		s *= "with $(length(midi.tracks)) tracks"
	end
	print(io, s)
end

MIDIFile() = MIDIFile(0,960,MIDITrack[])

"""
    qpm(midi)
Return the QPM (quarter notes per minute) where the given `MIDIFile` was exported at.
Returns 120 if not found.
"""
function qpm(t::MIDI.MIDIFile)
    # Find the one that corresponds to Set Tempo:
    # The event tttttt corresponds to the command
    # FF 51 03 tttttt Set Tempo (in microseconds per MIDI quarter-note)
    # See here (page 8):
    # http://www.cs.cmu.edu/~music/cmsip/readings/Standard-MIDI-file-format-updated.pdf
    for event in t.tracks[1].events
        if event isa SetTempoEvent
            return 6e7 / event.tempo
        end
    end

    # Default QPM if it is not present in the MIDI file.
    @warn """The Set Tempo event is not present in the given MIDI file.
    A default value of 120.0 quarter notes per minute is returned."""
    return 120.0
end

"""
    bpm(midi)
Return the BPM where the given `MIDIFile` was exported at.
Returns QPM if not found.
"""
function bpm(t::MIDI.MIDIFile)
    cc = -1

    # Find the one that corresponds to Time Signature:
    # FF 58 04 nn dd cc bb Time Signature
    # See here (page 8):
    # http://www.cs.cmu.edu/~music/cmsip/readings/Standard-MIDI-file-format-updated.pdf
    for event in t.tracks[1].events
        if event isa TimeSignatureEvent
            cc = event.clockticks
            break
        end
    end

    if cc == -1
        @warn """The Time Signature event is not present in the given MIDI file.
        A default value of 24 cc (clocks per metronome click) is used for calculating the BPM."""
        # Default cc if not found
        cc = 24
    end

    bpm = qpm(t) * 24 / cc
end

# Deprecated
"""
    BPM(midi)
Return the BPM where the given `MIDIFile` was exported at.
Returns 120 if not found.
"""
function BPM(t::MIDI.MIDIFile)
    @warn """This function is deprecated.
    It returns quarter notes per minute instead of beats per minute.
    Please use `bpm` for beats per minute and `qpm` for quarter notes per minute."""

    # Find the one that corresponds to Set-Time:
    # The event tttttt corresponds to the command
    # FF 51 03 tttttt Set Tempo (in microseconds per MIDI quarter-note)
    # See here (page 8):
    # http://www.cs.cmu.edu/~music/cmsip/readings/Standard-MIDI-file-format-updated.pdf
    for event in t.tracks[1].events
        if event isa SetTempoEvent
            return 6e7 / event.tempo
        end
    end

    # Default BPM if it is not present in the MIDI file.
    @warn """The Set Tempo event is not present in the given MIDI file.
    A default value of 120.0 quarter notes per minute is returned."""
    return 120.0
end

"""
    time_signature(midi)
Return the time signature of the given `MIDIFile`.
Returns 4/4 if it doesn't find a time signature.
"""
function time_signature(t::MIDI.MIDIFile)
    # Find the one that corresponds to Time Signature:
    # FF 58 04 nn dd cc bb Time Signature
    # See here (page 8):
    # http://www.cs.cmu.edu/~music/cmsip/readings/Standard-MIDI-file-format-updated.pdf
    for event in t.tracks[1].events
        if event isa TimeSignatureEvent
            ts = string(event.numerator) * "/" * string(event.denominator)
            return ts
        end
    end

    @warn """The Time Signature event is not present in the given MIDI file.
    A default value of 4/4 is returned."""

    # Default time signature if it is not present in the file
    return "4/4"
end

"""
    tempochanges(midi)
Return a vector of (position, tempo) tuples for all the tempo events in the given `MIDIFile`
where position is in absolute time (from the beginning of the file) in ticks
and tempo is in quarter notes per minute.
Returns [(0, 120.0)] if there are no tempo events.
"""
function tempochanges(midi::MIDIFile)
    # Stores (position, tempo) pairs
    # Calls qpm() to store the first tempo value
    # If there is no tempo event, qpm will warn and return 120.0
    tempo_changes = [(0, qpm(midi))]
    position = 0
    for event in midi.tracks[1].events
        position += event.dT
        if event isa SetTempoEvent
            qpm = 6e7 / event.tempo

            # Allow only one tempo change at the beginning
            if position == 0
                tempo_changes = [(0, qpm)]
            else
                push!(tempo_changes, (position, qpm))
            end
        end
    end

    tempo_changes
end

"""
    ms_per_tick(tpq, qpm)
    ms_per_tick(midi::MIDIFile)
Return how many miliseconds is one tick, based
on the quarter notes per minute `qpm` and ticks per quarter note `tpq`.
"""
ms_per_tick(midi::MIDI.MIDIFile, qpm = qpm(midi)) = ms_per_tick(midi.tpq, qpm)
ms_per_tick(tpq, qpm) = (1000*60)/(qpm*tpq)

getnotes(midi::MIDIFile, trackno = midi.format == 0 ? 1 : 2) =
getnotes(midi.tracks[trackno], midi.tpq)


text:
module child

f1() = 1
f2() = 2

end

ext: jl


text:
########### Defines multiplicaiton of SMC and SMC
function sigu(x::T) where {T<:AbstractFloat}
	 (zero(T)<=x) ? x^(1.0/MC_param.mu) : -abs(x)^(1.0/MC_param.mu)
 end
function xstar(y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
	return lambda.lo+(lambda.hi-lambda.lo)*(((nu.hi-y)/(nu.hi-nu.lo))
				 +sigu(-(nu.lo+nu.hi)/((MC_param.mu+1)*(nu.hi-nu.lo))))
end

function ystar(x::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
	return nu.lo+(nu.hi-nu.lo)*(((lambda.hi-x)/(lambda.hi-lambda.lo))
				 -sigu((lambda.lo+lambda.hi)/((MC_param.mu+1)*(lambda.hi-lambda.lo))))
end
@inline function dxA(x::T,y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
	      return (lambda.hi-lambda.lo)*(nu.hi-nu.lo)*
        abs((y-nu.lo)/(nu.hi-nu.lo)-(lambda.hi-x)/(lambda.hi-lambda.lo))^(MC_param.mu+1)
end
@inline function GxA(x::T,y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
	      return (x*(nu.lo+nu.hi)+y*(lambda.lo+lambda.hi)-
        (lambda.lo*nu.lo+lambda.hi*nu.hi)+dxA(x,y,lambda,nu))/2
end
@inline function dxB(x::T,y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
	      return (lambda.hi-lambda.lo)*(nu.hi-nu.lo)*max(zero(T),(((y-nu.lo)/(nu.hi-nu.lo))-(lambda.hi-x)/(lambda.hi-lambda.lo)))^(Int64(MC_param.mu+1))
end
@inline function GxB(x::T,y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
	      return x*nu.lo+y*lambda.lo-lambda.lo*nu.lo+dxB(x,y,lambda,nu)
end
@inline function gCxA(alpha::Interval{T},beta::Interval{T},lambda::Interval{T},
					  nu::Interval{T},x1::SMCg{N,T},x2::SMCg{N,T}) where {N,T<:AbstractFloat}
        term1::T,blank::Int64 = mid3(alpha.lo,alpha.hi,xstar(beta.lo,lambda,nu))
	    term2::T,blank = mid3(alpha.lo,alpha.hi,xstar(beta.hi,lambda,nu))
	    term3::T,blank = mid3(beta.lo,beta.hi,ystar(alpha.lo,lambda,nu))
	    term4::T,blank = mid3(beta.lo,beta.hi,ystar(alpha.hi,lambda,nu))
	    a::T,b::Int64 = findmin([GxA(term1,beta.lo,lambda,nu),GxA(term2,beta.hi,lambda,nu),
                           GxA(alpha.lo,term3,lambda,nu),GxA(alpha.hi,term4,lambda,nu)])
        		if (b == 1)
					#println("gCxA trace 1")
			  		grad::SVector{N,T} = max(zero(T),psi_mlt_Ax(term1,beta.lo,lambda,nu))*x1.cv_grad+
					   				 min(zero(T),psi_mlt_Ax(term1,beta.lo,lambda,nu))*x1.cc_grad+
					   	 			 max(zero(T),psi_mlt_Ay(term1,beta.lo,lambda,nu))*x2.cv_grad+
				 	   				 min(zero(T),psi_mlt_Ay(term1,beta.lo,lambda,nu))*x2.cc_grad

				    grad2::SVector{N,T} = max(zero(T),psi_mlt_Ax(term1,beta.lo,lambda,nu))*x1.cc_grad+
							min(zero(T),psi_mlt_Ax(term1,beta.lo,lambda,nu))*x1.cv_grad-
							max(zero(T),psi_mlt_Ay(term1,beta.lo,lambda,nu))*x2.cc_grad-
							min(zero(T),psi_mlt_Ay(term1,beta.lo,lambda,nu))*x2.cv_grad
				elseif (b == 2)
					#println("gCxA trace 2")
					grad = max(zero(T),psi_mlt_Ax(term2,beta.hi,lambda,nu))*x1.cv_grad+
						   min(zero(T),psi_mlt_Ax(term2,beta.hi,lambda,nu))*x1.cc_grad+
					 	   max(zero(T),psi_mlt_Ay(term2,beta.hi,lambda,nu))*x2.cv_grad+
				 		   min(zero(T),psi_mlt_Ay(term2,beta.hi,lambda,nu))*x2.cc_grad
				    grad2 = max(zero(T),psi_mlt_Ax(term2,beta.hi,lambda,nu))*x1.cc_grad+
							min(zero(T),psi_mlt_Ax(term2,beta.hi,lambda,nu))*x1.cv_grad-
							max(zero(T),psi_mlt_Ay(term2,beta.hi,lambda,nu))*x2.cc_grad-
							min(zero(T),psi_mlt_Ay(term2,beta.hi,lambda,nu))*x2.cv_grad
				elseif (b == 3)
				#	println("gCxA trace 3")
					grad = max(zero(T),psi_mlt_Ax(alpha.lo,term3,lambda,nu))*x1.cv_grad+
						   min(zero(T),psi_mlt_Ax(alpha.lo,term3,lambda,nu))*x1.cc_grad+
					 	   max(zero(T),psi_mlt_Ay(alpha.lo,term3,lambda,nu))*x2.cv_grad+
				 		   min(zero(T),psi_mlt_Ay(alpha.lo,term3,lambda,nu))*x2.cc_grad
				    grad2 = max(zero(T),psi_mlt_Ax(alpha.lo,term3,lambda,nu))*x1.cc_grad+
							min(zero(T),psi_mlt_Ax(alpha.lo,term3,lambda,nu))*x1.cv_grad-
							max(zero(T),psi_mlt_Ay(alpha.lo,term3,lambda,nu))*x2.cc_grad-
							min(zero(T),psi_mlt_Ay(alpha.lo,term3,lambda,nu))*x2.cv_grad
				else
				#	println("gCxA trace 4")
					grad = max(zero(T),psi_mlt_Ax(alpha.hi,term4,lambda,nu))*x1.cv_grad+
						   min(zero(T),psi_mlt_Ax(alpha.hi,term4,lambda,nu))*x1.cc_grad+
					 	   max(zero(T),psi_mlt_Ay(alpha.hi,term4,lambda,nu))*x2.cv_grad+
				 		   min(zero(T),psi_mlt_Ay(alpha.hi,term4,lambda,nu))*x2.cc_grad
				   grad2 = max(zero(T),psi_mlt_Ax(alpha.hi,term4,lambda,nu))*x1.cc_grad+
				        	min(zero(T),psi_mlt_Ax(alpha.hi,term4,lambda,nu))*x1.cv_grad-
							max(zero(T),psi_mlt_Ay(alpha.hi,term4,lambda,nu))*x2.cc_grad-
							min(zero(T),psi_mlt_Ay(alpha.hi,term4,lambda,nu))*x2.cv_grad
				end

        return a,grad,grad2
end

# include correction from Khan
function psi_mlt_Ax(x::T,y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
	term = [(y-nu.lo)/(nu.hi-nu.lo) (lambda.hi-x)/(lambda.hi-lambda.lo)]
	return 0.5*(nu.lo+nu.hi+(MC_param.mu+1)*(nu.hi-nu.lo)*(term[1]-term[2])*abs(term[1]-term[2])^(MC_param.mu-1))
end
# include correction from Khan
function psi_mlt_Ay(x::T,y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
	term = [(y-nu.lo)/(nu.hi-nu.lo) (lambda.hi-x)/(lambda.hi-lambda.lo)]
	return 0.5*(lambda.lo+lambda.hi+(MC_param.mu+1)*(lambda.hi-lambda.lo)*(term[1]-term[2])*abs(term[1]-term[2])^(MC_param.mu-1))
end
function psi_mlt_Bx(x::T,y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
	term = [(y - nu.lo)/(nu.hi-nu.lo)  (lambda.hi - x)/(lambda.hi-lambda.lo)]
	return nu.lo+(MC_param.mu+1)*(nu.hi-nu.lo)*max(0,term[1]-term[2])^MC_param.mu
end
function psi_mlt_By(x::T,y::T,lambda::Interval{T},nu::Interval{T}) where {T<:AbstractFloat}
	term = [(y - nu.lo)/(nu.hi-nu.lo)  (lambda.hi - x)/(lambda.hi-lambda.lo)]
	return lambda.lo+(MC_param.mu+1)*(lambda.hi-lambda.lo)*max(0,term[1]-term[2])^MC_param.mu
end

function multiply_MV(x1::SMCg{N,T},x2::SMCg{N,T}) where {N,T<:AbstractFloat}
	cv::T = zero(T)
	cc::T = zero(T)
	alpha0::Interval{T} = Interval(x1.cv,x1.cc)
	beta0::Interval{T} =  Interval(x2.cv,x2.cc)
	if (zero(T)<=x1.Intv.lo) && (zero(T)<=x2.Intv.lo)
		cv = GxB(x1.cv,x2.cv,x1.Intv,x2.Intv)
		cv_grad::SVector{N,T} = x1.cv_grad*psi_mlt_Bx(x1.cv,x2.cv,x1.Intv,x2.Intv) + x2.cv_grad*psi_mlt_By(x1.cv,x2.cv,x1.Intv,x2.Intv)
	elseif ((x1.Intv.hi<=zero(T))) && (x2.Intv.hi<=zero(T))
		cv = GxB(-x1.cc,-x2.cc,-x1.Intv,-x2.Intv)
		cv_grad = - x1.cc_grad*psi_mlt_Bx(-x1.cc,-x2.cc,-x1.Intv,-x2.Intv) - x2.cc_grad*psi_mlt_By(-x1.cc,-x2.cc,-x1.Intv,-x2.Intv)
	else
		cv,cv_grad,temp = gCxA(alpha0,beta0,x1.Intv,x2.Intv,x1,x2)
	end
	if ((x1.Intv.hi<=zero(T))) && ((zero(T))<=x2.Intv.lo)
		cc = -GxB(-x1.cc,x2.cv,-x1.Intv,x2.Intv)
		cc_grad::SVector{N,T} = x1.cc_grad*psi_mlt_Bx(-x1.cc,x2.cv,-x1.Intv,x2.Intv) - x2.cv_grad*psi_mlt_By(-x1.cc,x2.cv,-x1.Intv,x2.Intv)
	elseif ((zero(T))<=x1.Intv.lo) && (x2.Intv.hi<=zero(T))
		cc = -GxB(x1.cv,-x2.cc,x1.Intv,-x2.Intv)
		cc_grad = - x1.cv_grad*psi_mlt_Bx(x1.cv,-x2.cc,x1.Intv,-x2.Intv) + x2.cc_grad*psi_mlt_By(x1.cv,-x2.cc,x1.Intv,-x2.Intv)
	else
		cct::T,temp,cc_gradt::SVector{N,T} = gCxA(-alpha0,beta0,-x1.Intv,x2.Intv,x1,x2)
		cc = -cct
		cc_grad = cc_gradt
	end
	if (min(x1.Intv.lo,x2.Intv.lo)<zero(T)<max(x1.Intv.hi,x2.Intv.hi))
		lo_Intv_calc::T,blank = gCxA(x1.Intv,x2.Intv,x1.Intv,x2.Intv,x1,x2)
		hi_Intv_calct::T,blankt = gCxA(-x1.Intv,x2.Intv,-x1.Intv,x2.Intv,x1,x2)
		hi_Intv_calc::T = -hi_Intv_calct
		Intv_calc::Interval{T} = @interval(lo_Intv_calc,hi_Intv_calc)
	else
		Intv_calc = x1.Intv*x2.Intv
	end
	cnst = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
	return SMCg{N,T}(cc, cv, cc_grad, cv_grad, Intv_calc, cnst, x1.IntvBox,x1.xref)
end

function mul1_u1pos_u2pos(x1::SMCg{N,T},x2::SMCg{N,T},cnst::Bool) where {N,T<:AbstractFloat}
  Intv::Interval{T} = x1.Intv*x2.Intv
  cv1::T = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
  cv2::T = x2.Intv.lo*x1.cv + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
  if (cv1 > cv2)
    cv::T = cv1
    cv_grad::SVector{N,T} = x2.Intv.hi*x1.cv_grad
  else
    cv = cv2
    cv_grad = x2.Intv.lo*x1.cv_grad
  end

  cc1::T = x2.Intv.lo*x1.cc + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
  cc2::T = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
  if (cc1 < cc2)
    cc::T = cc1
    cc_grad::SVector{N,T} = x2.Intv.lo*x1.cc_grad
  else
    cc = cc2
    cc_grad = x2.Intv.hi*x1.cc_grad
  end
  return SMCg{N,T}(cc,cv,cc_grad,cv_grad,Intv,cnst,x1.IntvBox,x1.xref)
end
function mul1_u1pos_u2mix(x1::SMCg{N,T},x2::SMCg{N,T},cnst::Bool) where {N,T<:AbstractFloat}
  Intv::Interval{T} = x1.Intv*x2.Intv
  cv1::T = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
  cv2::T = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
  if (cv1 > cv2)
    cv::T = cv1
    cv_grad::SVector{N,T} = x2.Intv.hi*x1.cv_grad
  else
    cv = cv2
    cv_grad = x2.Intv.lo*x1.cc_grad
  end

  cc1::T = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
  cc2::T = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
  if (cc1 < cc2)
    cc::T = cc1
    cc_grad::SVector{N,T} = x2.Intv.lo*x1.cv_grad
  else
    cc = cc2
    cc_grad = x2.Intv.hi*x1.cc_grad
  end
  return SMCg{N,T}(cc,cv,cc_grad,cv_grad,Intv,cnst,x1.IntvBox,x1.xref)
end
function mul1_u1mix_u2mix(x1::SMCg{N,T},x2::SMCg{N,T},cnst::Bool) where {N,T<:AbstractFloat}
  Intv::Interval{T} = x1.Intv*x2.Intv
  cv1::T = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
  cv2::T = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.lo
  cv::T = cv1
  if (cv1 > cv2)
    cv = cv1
    cv_grad::SVector{N,T} = x2.Intv.hi*x1.cv_grad
  else
    cv = cv2
    cv_grad = x2.Intv.lo*x1.cc_grad
  end
  cc1::T = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
  cc2::T = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.hi
  cc::T = cc1
  if (cc1 < cc2)
    cc = cc1
    cc_grad::SVector{N,T} = x2.Intv.lo*x1.cv_grad
  else
    cc = cc2
    cc_grad = x2.Intv.hi*x1.cc_grad
  end
  return SMCg{N,T}(cc,cv,cc_grad,cv_grad,Intv,cnst,x1.IntvBox,x1.xref)
end
function mul2_u1pos_u2pos(x1::SMCg{N,T},x2::SMCg{N,T}) where {N,T<:AbstractFloat}
	Intv::Interval{T} = x1.Intv*x2.Intv
	cnst::Bool = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
	cv1::T = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
	cv2::T = x2.Intv.lo*x1.cv + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
	if (cv1 > cv2)
		cv::T = cv1
		cv_grad::SVector{N,T} = x2.Intv.hi*x1.cv_grad + x1.Intv.hi*x2.cv_grad
	else
		cv = cv2
		cv_grad = x2.Intv.lo*x1.cv_grad + x1.Intv.lo*x2.cv_grad
	end

	cc1::T = x2.Intv.lo*x1.cc + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
	cc2::T = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
	if (cc1 < cc2)
		cc::T = cc1
		cc_grad::SVector{N,T} = x2.Intv.lo*x1.cc_grad + x1.Intv.hi*x2.cc_grad
	else
		cc = cc2
		cc_grad = x2.Intv.hi*x1.cc_grad + x1.Intv.lo*x2.cc_grad
	end
	return SMCg{N,T}(cc,cv,cc_grad,cv_grad,Intv, cnst,x1.IntvBox,x1.xref)
end
function mul2_u1pos_u2mix(x1::SMCg{N,T},x2::SMCg{N,T},cnst::Bool) where {N,T<:AbstractFloat}
  Intv::Interval{T} = x1.Intv*x2.Intv
  cv1::T = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
  cv2::T = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
  if (cv1 > cv2)
    cv::T = cv1
    cv_grad::SVector{N,T} = x1.Intv.hi*x2.cv_grad
  else
    cv = cv2
    cv_grad = x1.Intv.lo*x2.cc_grad
  end

  cc1::T = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
  cc2::T = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
  if (cc1 < cc2)
    cc::T = cc1
    cc_grad::SVector{N,T} = x1.Intv.hi*x2.cc_grad
  else
    cc = cc2
    cc_grad = x1.Intv.lo*x2.cc_grad
  end
  return SMCg{N,T}(cc,cv,cc_grad,cv_grad,Intv,cnst,x1.IntvBox,x1.xref)
end
function mul2_u1mix_u2mix(x1::SMCg{N,T},x2::SMCg{N,T}) where {N,T<:AbstractFloat}
	Intv::Interval{T} = x1.Intv*x2.Intv
  	cnst::Bool = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
	cv1::T = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
	cv2::T = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.lo
	if (cv1 > cv2)
		cv::T = cv1
		cv_grad::SVector{N,T} = x2.Intv.hi*x1.cv_grad + x1.Intv.hi*x2.cv_grad
	else
		cv = cv2
		cv_grad = x2.Intv.lo*x1.cc_grad + x1.Intv.lo*x2.cc_grad
	end

	cc1::T = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
	cc2::T = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.hi
	if (cc1 < cc2)
		cc::T = cc1
		cc_grad::SVector{N,T} = x2.Intv.lo*x1.cv_grad + x1.Intv.hi*x2.cc_grad
	else
		cc = cc2
		cc_grad = x2.Intv.hi*x1.cc_grad + x1.Intv.lo*x2.cv_grad
	end

	return SMCg{N,T}(cc,cv,cc_grad,cv_grad,Intv,cnst,x1.IntvBox,x1.xref)
end
function mul3_u1pos_u2mix(x1::SMCg{N,T},x2::SMCg{N,T}) where {N,T<:AbstractFloat}
	Intv::Interval{T} = x1.Intv*x2.Intv
    cnst::Bool = x2.cnst ? x1.cnst : (x1.cnst ? x2.cnst : x1.cnst || x2.cnst)
	cv1::T = x2.Intv.hi*x1.cv + x1.Intv.hi*x2.cv - x1.Intv.hi*x2.Intv.hi
	cv2::T = x2.Intv.lo*x1.cc + x1.Intv.lo*x2.cv - x1.Intv.lo*x2.Intv.lo
	if (cv1 > cv2)
		cv::T = cv1
		cv_grad::SVector{N,T} = x2.Intv.hi*x1.cv_grad + x1.Intv.hi*x2.cv_grad
	else
		cv = cv2
		cv_grad = x2.Intv.lo*x1.cc_grad + x1.Intv.lo*x2.cv_grad
	end

	cc1::T = x2.Intv.lo*x1.cv + x1.Intv.hi*x2.cc - x1.Intv.hi*x2.Intv.lo
	cc2::T = x2.Intv.hi*x1.cc + x1.Intv.lo*x2.cc - x1.Intv.lo*x2.Intv.hi
	if (cc1 < cc2)
		cc::T = cc1
		cc_grad::SVector{N,T} = x2.Intv.lo*x1.cv_grad + x1.Intv.hi*x2.cc_grad
	else
		cc = cc2
		cc_grad = x2.Intv.hi*x1.cc_grad + x1.Intv.lo*x2.cc_grad
	end

	return SMCg{N,T}(cc,cv,cc_grad,cv_grad,Intv,cnst,x1.IntvBox,x1.xref)
end

function mul_MV_ns1cv(x1::T,x2::T,MC1::SMCg{N,T},MC2::SMCg{N,T}) where {N,T<:AbstractFloat}
	return MC2.Intv.hi*x1+MC1.Intv.hi*x2-MC2.Intv.hi*MC1.Intv.hi
end
function mul_MV_ns2cv(x1::T,x2::T,MC1::SMCg{N,T},MC2::SMCg{N,T}) where {N,T<:AbstractFloat}
	return MC2.Intv.lo*x1+MC1.Intv.lo*x2-MC2.Intv.lo*MC1.Intv.lo
end
function mul_MV_ns3cv(x1::T,x2::T,MC1::SMCg{N,T},MC2::SMCg{N,T}) where {N,T<:AbstractFloat}
	return max(mul_MV_ns1cv(x1,x2,MC1,MC2),mul_MV_ns2cv(x1,x2,MC1,MC2))
end
function mul_MV_ns1cc(x1::T,x2::T,MC1::SMCg{N,T},MC2::SMCg{N,T}) where {N,T<:AbstractFloat}
	return MC2.Intv.lo*x1+MC1.Intv.hi*x2-MC2.Intv.lo*MC1.Intv.hi
end
function mul_MV_ns2cc(x1::T,x2::T,MC1::SMCg{N,T},MC2::SMCg{N,T}) where {N,T<:AbstractFloat}
	return MC2.Intv.hi*x1+MC1.Intv.lo*x2-MC2.Intv.hi*MC1.Intv.lo
end
function mul_MV_ns3cc(x1::T,x2::T,MC1::SMCg{N,T},MC2::SMCg{N,T}) where {N,T<:AbstractFloat}
	return min(mul_MV_ns1cc(x1,x2,MC1,MC2),mul_MV_ns2cc(x1,x2,MC1,MC2))
end
function multiply_MV_NS(x1::SMCg{N,T},x2::SMCg{N,T},ngrad::Int64,cnst::Bool) where {N,T<:AbstractFloat}

 k::T = diam(x2.Intv)/diam(x1.Intv)
 z::T = (x1.Intv.hi*x2.Intv.hi - x1.Intv.lo*x2.Intv.lo)/diam(x1.Intv)
 x1vta::T,blank::Int64 = mid3(x1.cv,x1.cc,(x2.cv-z)/k)
 x1vtb::T,blank = mid3(x1.cv,x1.cc,(x2.cc-z)/k)
 x2vta::T,blank = mid3(x2.cv,x2.cc, k*x1.cv+z)
 x2vtb::T,blank = mid3(x2.cv,x2.cc, k*x1.cc+z)
 x1vt = [x1.cv, x1.cc, x1vta, x1vtb, x1.cv, x1.cc]
 x2vt = [x2vta, x2vtb, x2.cv, x2.cc, x2.cv, x2.cc]
 vt  = [mul_MV_ns3cv(x1vt[1],x2vt[1],x1,x2), mul_MV_ns3cv(x1vt[2],x2vt[2],x1,x2),
        mul_MV_ns3cv(x1vt[3],x2vt[3],x1,x2), mul_MV_ns3cv(x1vt[4],x2vt[4],x1,x2),
				mul_MV_ns3cv(x1vt[5],x2vt[5],x1,x2), mul_MV_ns3cv(x1vt[6],x2vt[6],x1,x2)]
 cv::T,cvind::Int64 = findmax(vt)

 if (ngrad>0)
	if isequal(mul_MV_ns1cv(x1vt[cvind],x2vt[cvind],x1,x2),
						 mul_MV_ns2cv(x1vt[cvind],x2vt[cvind],x1,x2),
						 MC_param.mv_tol,MC_param.mv_tol)
 		alph = [zero(T),one(T)]

 		MC1thin::Bool = isequal(x1.cv,x1.cc,MC_param.mv_tol,MC_param.mv_tol)
 		if ((~MC1thin) && (x1vt[cvind] > x1.cv))
	 		if (~isequal(x1vt[cvind],x1.cv,MC_param.mv_tol,MC_param.mv_tol))
		 		alph[2] = min(alph[2],-x2.Intv.lo/diam(x2.Intv))
	 		end
 		end
 		if ((~MC1thin) && (x1vt[cvind] < x1.cc))
			if (~isequal(x1vt[cvind],x1.cc,MC_param.mv_tol,MC_param.mv_tol))
				alph[1] = max(alph[1],-x2.Intv.lo/diam(x2.Intv))
			end
 		end

 		MC2thin::Bool = isequal(x2.cv,x2.cc,MC_param.mv_tol,MC_param.mv_tol)
 		if ((~MC2thin) && (x2vt[cvind] > x2.cv))
			if (~isequal(x2vt[cvind],x2.cv,MC_param.mv_tol,MC_param.mv_tol))
				alph[2] = min(alph[2],-x1.Intv.lo/diam(x1.Intv))
			end
 		end
 		if ((~MC2thin) && (x2vt[cvind] < x2.cc))
 			if (~isequal(x2vt[cvind],x2.cc,MC_param.mv_tol,MC_param.mv_tol))
	 			alph[1] = max(alph[1],-x1.Intv.lo/diam(x1.Intv))
 			end
 		end

 		alphthin::Bool = isequal(alph[1],alph[2],MC_param.mv_tol,MC_param.mv_tol)
 		if (~alphthin && (alph[1]>alph[2]))
	 		error("Multivariant mult error alphaL = alphaU")
 		end
 		myalph::T = (alph[1]+alph[2])/2
 	elseif (mul_MV_ns1cv(x1vt[cvind],x2vt[cvind],x1,x2) >
			    mul_MV_ns2cv(x1vt[cvind],x2vt[cvind],x1,x2))
		myalph = one(T)
	else
		myalph = zero(T)
	end
	sigma_cv1::T = x2.Intv.lo + myalph*diam(x2.Intv)
	sigma_cv2::T = x1.Intv.lo + myalph*diam(x1.Intv)
	if (x1.cnst)
		term1::SVector{N,T} = @SVector zeros(T,N)
	elseif (sigma_cv1>=zero(T))
		term1 = x1.cv_grad
	else
		term1 = x1.cc_grad
	end
	if (x2.cnst)
		term2::SVector{N,T} = @SVector zeros(T,N)
	elseif (sigma_cv1>=zero(T))
		term2 = x2.cv_grad
	else
		term2 = x2.cc_grad
	end
	cv_grad::SVector{N,T} = term1*sigma_cv1 + term2*sigma_cv2
 end

 z = (x1.Intv.hi*x2.Intv.lo - x1.Intv.lo*x2.Intv.hi)/diam(x1.Intv)
 x1vta,blank = mid3(x1.cv,x1.cc,(x2.cv-z)/k)
 x1vtb,blank = mid3(x1.cv,x1.cc,(x2.cc-z)/k)
 x2vta,blank = mid3(x2.cv,x2.cc, k*x1.cv+z)
 x2vtb,blank = mid3(x2.cv,x2.cc, k*x1.cc+z)
 x1vt = [x1.cv, x1.cc, x1vta, x1vtb, x1.cv, x1.cc]
 x2vt = [x2vta, x2vtb, x2.cv, x2.cc, x2.cc, x2.cv]
 vt  = [mul_MV_ns3cc(x1vt[1],x2vt[1],x1,x2), mul_MV_ns3cc(x1vt[2],x2vt[2],x1,x2),
				mul_MV_ns3cc(x1vt[3],x2vt[3],x1,x2), mul_MV_ns3cc(x1vt[4],x2vt[4],x1,x2),
			  mul_MV_ns3cc(x1vt[5],x2vt[5],x1,x2), mul_MV_ns3cc(x1vt[6],x2vt[6],x1,x2)]
 cc::T,ccind::Int64 = findmax(vt)

 if (ngrad>0)
 	if isequal(mul_MV_ns1cc(x1vt[cvind],x2vt[cvind],x1,x2),
						 mul_MV_ns2cc(x1vt[cvind],x2vt[cvind],x1,x2),
						 MC_param.mv_tol,MC_param.mv_tol)
		 alph = [zero(T),one(T)]

		 MC1thin = isequal(x1.cv,x1.cc,MC_param.mv_tol,MC_param.mv_tol)
		 if ((~MC1thin) && (x1vt[cvind] > x1.cv))
		 	if (~isequal(x1vt[cvind],x1.cv,MC_param.mv_tol,MC_param.mv_tol))
			 	alph[1] = max(alph[1],-x2.Intv.lo/diam(x2.Intv))
		 	end
		 end
		 if ((~MC1thin) && (x1vt[cvind] < x1.cc))
		 	if (~isequal(x1vt[cvind],x1.cc,MC_param.mv_tol,MC_param.mv_tol))
			 	alph[2] = min(alph[2],-x2.Intv.lo/diam(x2.Intv))
		 	end
		 end

		 MC2thin = isequal(x2.cv,x2.cc,MC_param.mv_tol,MC_param.mv_tol)
		 if ((~MC2thin) && (x2vt[cvind] > x2.cv))
		 	if (~isequal(x2vt[cvind],x2.cv,MC_param.mv_tol,MC_param.mv_tol))
			 	alph[2] = min(alph[2],x1.Intv.hi/diam(x1.Intv))
		 	end
		 end
		 if ((~MC2thin) && (x2vt[cvind] < x2.cc))
			if (~isequal(x2vt[cvind],x2.cc,MC_param.mv_tol,MC_param.mv_tol))
			 	alph[1] = max(alph[1],x1.Intv.hi/diam(x1.Intv))
			end
		 end

		 alphthin = isequal(alph[1],alph[2],MC_param.mv_tol,MC_param.mv_tol)
		 if (~alphthin && (alph[1]>alph[2]))
		 	error("Multivariant mult error alphaL = alphaU")
		 end
		 myalph = (alph[1]+alph[2])/2
	elseif (mul_MV_ns1cv(x1vt[cvind],x2vt[cvind],x1,x2) >
				  mul_MV_ns2cv(x1vt[cvind],x2vt[cvind],x1,x2))
	   myalph = one(T)
  	else
		myalph = zero(T)
  	end
	sigma_cc1::T = x2.Intv.lo + myalph*diam(x2.Intv)
	sigma_cc2::T = x1.Intv.hi - myalph*diam(x1.Intv)
	if (x1.cnst)
		term1 = @SVector zeros(T,N)
	elseif (sigma_cc1>=zero(T))
		term1 = x1.cc_grad
	else
		term1 = x1.cv_grad
	end
	if (x2.cnst)
		term2 =  @SVector zeros(T,N)
	elseif (sigma_cc1>=zero(T))
		term2 = x2.cc_grad
	else
		term2 = x2.cv_grad
	end
	cc_grad = term1*sigma_cc1 + term2*sigma_cc2
 end
 return SMCg{N,T}(cc,cv,cc_grad,cv_grad,x1.Intv*x2.Intv,cnst,x1.IntvBox,x1.xref)
end

function multiply_STD_NS(x1::SMCg{N,T},x2::SMCg{N,T}) where {N,T<:AbstractFloat}
	if (x2.Intv.lo >= zero(T))
    	if (x2.cnst)
      		return mul1_u1pos_u2pos(x1,x2,x1.cnst) # cnst to do
    	elseif (x1.cnst)
      		return mul1_u1pos_u2pos(x2,x1,x2.cnst) # cnst to do
    	else
      		return mul2_u1pos_u2pos(x1,x2) # DONE
    	end
	elseif (x2.Intv.hi <= zero(T))
		return -(x1*(-x2))
	else
    	if (x2.cnst)
      		return mul1_u1pos_u2mix(x1,x2,x1.cnst) # cnst to do
    	elseif (x1.cnst)
      		return mul2_u1pos_u2mix(x1,x2,x2.cnst) # cnst to do
    	else
	  		return mul3_u1pos_u2mix(x1,x2) # cnst to do
    	end
	end
end

function STD_NS_ALT(x::SMCg{N,T},y::SMCg{N,T}) where {N,T<:AbstractFloat}
	alpha1::T = min( y.Intv.lo*x.cv,  y.Intv.lo*x.cc )
	alpha2::T = min( x.Intv.lo*y.cv,  x.Intv.lo*y.cc )
	beta1::T  = min( y.Intv.hi*x.cv,  y.Intv.hi*x.cc )
	beta2::T  = min( x.Intv.hi*y.cv,  x.Intv.hi*y.cc )
	gamma1::T = max( y.Intv.lo*x.cv,  y.Intv.lo*x.cc )
	gamma2::T = max( x.Intv.hi*y.cv,  x.Intv.hi*y.cc )
	delta1::T = max( y.Intv.hi*x.cv,  y.Intv.hi*x.cc )
	delta2::T = max( x.Intv.lo*y.cv,  x.Intv.lo*y.cc )

	cv1::T = alpha1 + alpha2 - x.Intv.lo*y.Intv.lo
	cv2::T = beta1  + beta2  - x.Intv.hi*y.Intv.hi
	cc1::T = gamma1 + gamma2 - x.Intv.hi*y.Intv.lo
	cc2::T = delta1 + delta2 - x.Intv.lo*y.Intv.hi

	s_alpha1::SVector{N,T} = (y.Intv.lo >= zero(T)) ? y.Intv.lo*x.cv_grad : y.Intv.lo*x.cc_grad
	s_alpha2::SVector{N,T} = (x.Intv.lo >= zero(T)) ? x.Intv.lo*x.cv_grad : x.Intv.lo*x.cc_grad
	s_beta1::SVector{N,T}  = (y.Intv.hi >= zero(T)) ? y.Intv.hi*x.cv_grad : y.Intv.hi*x.cc_grad
	s_beta2::SVector{N,T}  = (x.Intv.hi >= zero(T)) ? x.Intv.hi*x.cv_grad : x.Intv.hi*x.cc_grad
	s_gamma1::SVector{N,T} = (y.Intv.lo >= zero(T)) ? y.Intv.lo*x.cc_grad : y.Intv.lo*x.cv_grad
	s_gamma2::SVector{N,T} = (x.Intv.hi >= zero(T)) ? x.Intv.hi*x.cc_grad : x.Intv.hi*x.cv_grad
	s_delta1::SVector{N,T} = (y.Intv.hi >= zero(T)) ? y.Intv.hi*x.cc_grad : y.Intv.hi*x.cv_grad
	s_delta2::SVector{N,T} = (x.Intv.lo >= zero(T)) ? x.Intv.lo*x.cc_grad : x.Intv.lo*x.cv_grad

	if (cv1 >= cv2)
	    cv::T = cv1
	    cv_grad::SVector{N,T} = s_alpha1 + s_alpha2
	else
	    cv = cv2
	    cv_grad = s_beta1 + s_beta2
	end

	if (cc1 <= cc2)
	    cc::T = cc1
	    cc_grad::SVector{N,T} = s_gamma1 + s_gamma2
	else
	    cc = cc2
	    cc_grad = s_delta1 + s_delta2
	end
	Intv::Interval{T} = x.Intv*y.Intv
	return SMCg{N,T}(cc,cv,cc_grad,cv_grad,Intv,(x.cnst && y.cnst),x.IntvBox,x.xref)
end

@inline function *(x1::SMCg{N,T},x2::SMCg{N,T}) where {N,T<:AbstractFloat}
	if x1 == x2
		println("sqr trace")
		return sqr(x1)
	end

	degen1::Bool = ((x1.Intv.hi - x1.Intv.lo) == zero(T))
	degen2::Bool = ((x2.Intv.hi - x2.Intv.lo) == zero(T))

	if (MC_param.mu >= 1 && ~(degen1||degen2))
		return multiply_MV(x1,x2)
	elseif (MC_param.multivar_refine && ~(degen1||degen2))
		println("NS MV mult trace 1")
		if (x2.cnst)
			cnst = x1.cnst
		elseif (x1.cnst)
			cnst = x2.cnst
		elseif (length(x1.cc_grad) != length(x2.cc_grad))
			error("Unequal subgradients")
		else
			cnst = (x1.cnst||x2.cnst)
		end
		return multiply_MV_NS(x1,x2,N,cnst) # DONE (minus gradients & case handling?)
	elseif (x1.Intv.lo >= zero(T))
		return multiply_STD_NS(x1,x2)
		#return STD_NS_ALT(x1,x2)
	elseif (x1.Intv.hi <= zero(T))
		if (x2.Intv.lo >= zero(T))
			return -((-x1)*x2)
		elseif (x2.Intv.hi <= zero(T))
			return (-x1)*(-x2)
		else
			return -(x2*(-x1))
		end
	elseif (x2.Intv.lo >= zero(T))
		return x2*x1
	elseif (x2.Intv.hi <= zero(T))
		return -((-x2)*x1)
	else
    	if (x2.cnst)
			#println("NS mult trace 5")
	  		return STD_NS_ALT(x1,x2)
      		#return mul1_u1mix_u2mix(x1,x2,x1.cnst)
    	elseif (x1.cnst)
	  		#println("NS mult trace 6")
	  		return STD_NS_ALT(x1,x2)
      		#return mul1_u1mix_u2mix(x2,x1,x2.cnst)
    	else
	  		#return STD_NS_ALT(x1,x2)
	  		return mul2_u1mix_u2mix(x1,x2)
    	end
	end
end


text:
# Util functions
include("./Classification.jl")

"""
    report_stats(y_pred, y_true; classwise_stats=true, avg_type="macro", sample_weights=nothing)

A utility function that prints the statistics summary of the model based on provided `y_pred` and `y_true`.

# Arguments:
 - `y_pred`: predicted values
 - `y_true`: ground truth values on the basis of which predicted values are to be assessed.
 - `classwise_stats=true`: if set `true`, prints classwise stats along with global stats.
 - `avg_type="macro"`: Type of average to be used while calculating precision of multiclass models. Can take values as `macro`, `micro` and `weighted`. Default set to `macro`.
 - `sample_weights`: Class weights to be provided when `avg_type` is set to `weighted`. Useful in case of imbalanced classes.

"""
function report_stats(y_pred, y_true; classwise_stats=true, avg_type="macro", sample_weights=nothing)
    g_stats = global_stats(y_pred, y_true, avg_type = avg_type)
    print("++++++++++++++++++++ Confusion Matrix +++++++++++++++++++\n")
    cm = g_stats["Confusion_Matrix"]
    for i in 1:size(cm,1)
        print(" ")
        for j in 1:size(cm,2)
            print(cm[i,j]," ")
        end
        print("\n")
    end
    print("\n")
    print("++++++++++++++++++++ Global Statistics ++++++++++++++++++++\n")
    print("  Accuracy: ", g_stats["Accuracy"],"\n")
    print("  Precision: ", g_stats["Precision"],"\n")
    print("  Recall: ", g_stats["Recall"],"\n")
    print("  F1_score: ", g_stats["F1_score"],"\n")
    print("  Specificity: ", g_stats["Specificity"],"\n")
    print("  False alarm rate: ", g_stats["False_alarm_rate"],"\n")
    print("\n")
    if (classwise_stats == true)
        c_stats = classwise_stats(y_pred, y_true)
        stats = DataFrame(Class=[], Confusion_Matrix=[], Accuracy=[], Precision=[], Recall=[], F1_score=[], Specificity=[])
        for i in 1: length(c_stats)
            push!(stats, Dict(:Class=>i,c_stats[i]...))
        end
        print(stats)
    end
end

# TODO report_stats_neat() - a function to print the Statistics summary more neatly

# TODO stats_to_csv() - function to save stats to a csv file.


text:
include("Hertz.jl")
include("Granular_Module.jl")
include("Fast_Hertz.jl")
#include("Multi_HZ.jl")
#include("calculate_Pressure.jl")
#include("Min_Energy.jl")
#include("Deform_Config.jl")

using JuLIP
#using Minimise_Energy
#using Deform
#using Pressure
using Hertz
using NeighbourLists
using Fast_Hertz
#using MultiHZ
# ===============================
#            new test
# ===============================

Grain_Num = 1024
α = 2.0
ϕ = 0.842

Jamming_Info = readdlm("jamming_config.txt")

x = Jamming_Info[:,1]'
y = Jamming_Info[:,2]'
z = zeros(Grain_Num)'
X = [x; y; z]
X = vecs(X)

grain_type = Jamming_Info[:,4] + 1
r = Jamming_Info[:,3]

Config = gran.setup_cell(Grain_Num, ϕ, α)
Config.X .= X
Config.M = r
Config.Z = grain_type

ϕ = π*sum(r.^2)
set_data!(Config, :Phi, ϕ)
set_data!(Config, :Rad, r)
Config1 = deepcopy(Config)

Num_Cal = 1000
# get the calculator
α = get_data(Config, :Alpha)
k = 1.0
r = Config.M
rmax = maximum(r)
rmin = minimum(r)
r0 = 2*rmax
r1 = 2*rmin

V = Hertz.HHZPotential(k, α, r0; rcutfact = 2.5)
set_calculator!(Config, V)

α1 = [α α; α α]
k1 = [k k; k k]
z = [1, 2]
#r1 = [2.0*rmin rmin + rmax; rmin + rmax 2.0*rmax]
V1 = Fast_Hertz.Fast_Hertz_Potential(get_data(Config, :Length), r0, r1; rcutfact = 2.5)
set_calculator!(Config1, V1)
E = 0.0
E1 = 0.0
tic()
 @simd for i = 1:Num_Cal #in range(1, Num_Cal)
    @inbounds E = Hertz.energy(V, Config)
      end
toc()

 tic()
  @simd for i = 1:Num_Cal #in range(1, Num_Cal)
     @inbounds E1 = Fast_Hertz.energy(V1, Config1)
       end
 toc()

# tic()
# for i = 1:Num_Cal
#     nlist = neighbourlist(Config, cutoff(V))
# end
# @show time2 = toc()

# tic()
# for i in range(1, Num_Cal)
#     E1 = energy(V1, Config1)
# end
# time2 = toc()


text:
using LinearAlgebra, Test
using CUDA

import LinearAlgebra: BlasInt

m = 256
n = 512

# our wrappers set up CUSOLVERMG for use with all available devices.
# on CI, if we only have a single device, set up multiple devices
# so that we properly cover the multigpu code paths.
if ndevices() == 1
    CUSOLVER.devices!([device(), device()])
end

@testset "mg_syevd!" begin
    @testset "element type $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
        A = rand(elty, m, m)
        A += A'
        hW = eigvals(Hermitian(A))
        hV = eigvecs(Hermitian(A))
        W, A = CUSOLVER.mg_syevd!('V','L',A)
        # compare
        @test W ≈ hW
        @test A*diagm(0=>W)*A' ≈ hV*diagm(0=>hW)*hV'

        A = rand(elty, m, m)
        A += A'
        hW = eigvals(Hermitian(A))
        W = CUSOLVER.mg_syevd!('N','L',A)
        # compare
        @test W ≈ hW

    end
end # elty

if CUSOLVER.version() >= v"10.4"
    @testset "mg_potrf!" begin
        @testset "element type $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
            A = rand(elty, m, m)
            A = A*A'+I # posdef
            hA = copy(A)
            A = CUSOLVER.mg_potrf!('L',A)
            LinearAlgebra.LAPACK.potrf!('L', hA)
            # compare
            @test A ≈ hA
        end
    end # elty

    @testset "mg_potrf and mg_potri!" begin
        #@testset "element type $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
        @testset "element type $elty" for elty in [Float64, ComplexF64]
            A = rand(elty, m, m)
            A = A*A'+I # posdef
            hA = copy(A)
            LinearAlgebra.LAPACK.potrf!('L', hA)
            A = CUSOLVER.mg_potrf!('L',A)
            LinearAlgebra.LAPACK.potri!('L', hA)
            A = CUSOLVER.mg_potri!('L',A)
            # compare
            @test tril(A) ≈ tril(hA)
        end
    end # elty

    @testset "mg_potrf and mg_potrs!" begin
        #@testset "element type $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
        @testset "element type $elty" for elty in [Float64, ComplexF64]
            A = rand(elty, m, m)
            B = rand(elty, m, m)
            A = A*A'+I # posdef
            hA = copy(A)
            hB = copy(B)
            LinearAlgebra.LAPACK.potrf!('L', hA)
            LinearAlgebra.LAPACK.potrs!('L', hA, hB)
            A = CUSOLVER.mg_potrf!('L', A)
            B = CUSOLVER.mg_potrs!('L', A, B)
            # compare
            tol    = real(elty) == Float32 ? 1e-4 : 1e-6
            @test A ≈ hA
            @test B ≈ hB rtol=tol
        end
    end # elty
end

if CUSOLVER.version() >= v"10.3"
    @testset "mg_getrf!" begin
        @testset "element type $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
            A      = rand(elty,m,m)
            h_A    = copy(A)
            A,ipiv = CUSOLVER.mg_getrf!(A)
            alu    = LinearAlgebra.LU(A, convert(Vector{BlasInt},ipiv), zero(BlasInt))
            @test h_A ≈ Array(alu)
        end
    end

    @testset "mg_getrs!" begin
        @testset "element type $elty" for elty in [Float32, Float64, ComplexF32, ComplexF64]
            A      = rand(elty,m,m)
            h_A    = copy(A)
            pivot = VERSION >= v"1.7-" ? NoPivot() : Val(false)
            alu    = lu(A, pivot)
            B      = rand(elty, m, div(m,2))
            h_B    = copy(B)
            tol    = real(elty) == Float32 ? 1e-1 : 1e-6
            B      = CUSOLVER.mg_getrs!('N', alu.factors, alu.ipiv, B)
            @test B ≈ h_A\h_B  rtol=tol
        end
    end
end


text:
# modified from https://github.com/chkwon/TrafficAssignment.jl

using CSV, DataFrames, LightGraphs, Optim, BinDeps;

type TA_Data
    network_name::String

    number_of_zones::Int
    number_of_nodes::Int
#     first_thru_node::Int
    number_of_links::Int

    start_node::Array{Int,1}
    end_node::Array{Int,1}
    capacity::Array{Float64,1}
    link_length::Array{Float64,1}
    free_flow_time::Array{Float64,1}
    B::Array{Float64,1}
    power::Array{Float64,1}
    speed_limit::Array{Float64,1}
    toll::Array{Float64,1}
    link_type::Array{Int64,1}

    total_od_flow::Float64

    travel_demand::Array{Float64,2}
    od_pairs::Array{Tuple{Int64,Int64},1}

    toll_factor::Float64
    distance_factor::Float64

    best_objective::Float64

    travel_time::Array{Float64,1}
end

net_filename = ARGS[1]
net_metadata_filename = ARGS[2]

trip_metadata_filename = ARGS[3]
trip_filename = ARGS[4]

output_filename = ARGS[5]

function load_ta_network(net_name, net_metadata_filename, net_filename, trip_metadata_filename, trip_filename;
    best_objective=-1.0, toll_factor=0.0, distance_factor=0.0, travel_time_initialization=nothing, fftime=:fftime)

    @assert ispath(net_metadata_filename)
    @assert ispath(net_filename)
    @assert ispath(trip_filename)

    ##################################################
    # Network Data
    ##################################################

    # read in metadata
    dfnm = CSV.read(net_metadata_filename; header=["Name", "Count"], allowmissing=:none)

    number_of_zones = dfnm[1,2]
    number_of_nodes = dfnm[2,2]
    #     first_thru_node = dfnm[3,2]
    number_of_links = dfnm[3,2]

    @assert number_of_links > 0

    # read in net
    dfn = CSV.read(net_filename; allowmissing=:none)

    start_node = dfn[:tail]
    end_node = dfn[:head]
    capacity = dfn[:capacity]
    link_length = dfn[:length]
    free_flow_time = dfn[fftime]
    B = dfn[:b]
    power = dfn[:power]
    speed_limit = dfn[:speedlimit]
    toll = dfn[:toll]
    link_type = dfn[:type]
    if travel_time_initialization != nothing
        travel_time = dfn[travel_time_initialization]
    else
        travel_time = zeros(number_of_links)
    end

    ##################################################
    # Trip Table
    ##################################################

    dftm = CSV.read(trip_metadata_filename; header=["Name", "Count"], allowmissing=:none)

    number_of_zones_trip = dftm[1,2]
    total_od_flow = dftm[2,2]

    @assert number_of_zones_trip == number_of_zones # Check if number_of_zone is same in both txt files
    @assert total_od_flow > 0

    travel_demand = zeros(number_of_zones, number_of_zones)
    od_pairs = Array{Tuple{Int64, Int64}}(0)

    dft = CSV.read(trip_filename; allowmissing=:none)

    ## old: columnar format
    # for row in eachrow(dft)
    #     push!(od_pairs, (row[1], row[2]))
    #     travel_demand[row[1], row[2]] = row[3]
    # end

    # new: od table format
    sort!(dft,1)
    delete!(dft,1)
    for i = 1:size(dft,2)
        for j = 1:size(dft,1)
            travel_demand[i,j] = dft[i,j]
        end
    end

    # Preparing data to return
    ta_data = TA_Data(
        net_name,
        number_of_zones,
        number_of_nodes,
    #         first_thru_node,
        number_of_links,
        start_node,
        end_node,
        capacity,
        link_length,
        free_flow_time,
        B,
        power,
        speed_limit,
        toll,
        link_type,
        total_od_flow,
        travel_demand,
        od_pairs,
        toll_factor,
        distance_factor,
        best_objective,
        travel_time
    )

    return ta_data

end # end of load_network function

dat = load_ta_network("moffet", net_metadata_filename, net_filename, trip_metadata_filename, trip_filename;
    travel_time_initialization=:am, fftime=:base)

function TA_dijkstra_shortest_paths(graph, travel_time, origin, start_node, end_node)
    no_node = nv(graph)
    no_arc = ne(graph)

    distmx = Inf*ones(no_node, no_node)
    for i in 1:no_arc
        distmx[start_node[i], end_node[i]] = travel_time[i]
    end

    state = dijkstra_shortest_paths(graph, origin, distmx)
    return state
end

function create_graph(start_node, end_node, number_of_nodes)
    @assert Base.length(start_node)==Base.length(end_node)

    no_arc = Base.length(start_node)

    graph = DiGraph(number_of_nodes)
    for i=1:no_arc
        add_edge!(graph, start_node[i], end_node[i])
    end
    return graph
end

function get_vector(state, origin, destination, link_dic)
    current = destination
    parent = -1
    x = zeros(Int, maximum(link_dic))
#     x = zeros(Int, Base.length(link_dic))

    while parent != origin && origin != destination && current != 0
        parent = state.parents[current]

        # println("origin=$origin, destination=$destination, parent=$parent, current=$current")

        if parent != 0
            link_idx = link_dic[parent,current]
            if link_idx != 0
                x[link_idx] = 1
            end
#             link_idx = get(link_dic, (parent,current), 0)
#             if link_idx != 0
#                 x[link_idx] = 1
#             end
        end

        current = parent
    end

    return x
end

function add_demand_vector!(x, demand, state, origin, destination, link_dic)
    current = destination
    parent = -1

    while parent != origin && origin != destination && current != 0
        parent = state.parents[current]

        if parent != 0
            link_idx = link_dic[parent,current]
            if link_idx != 0
                x[link_idx] += demand
            end
#           link_idx = get(link_dic, (parent,current), 0)
#           if link_idx != 0
#               x[link_idx] += demand
#           end
        end

        current = parent
    end
end

# Frank-Wolfe Methods
# CFW and BFW in Mitradijieva and Lindberg (2013)

# required packages: Graphs, Optim

# include("misc.jl")

function ta_frank_wolfe(ta_data; method=:bfw, min_iter_no=5, max_iter_no=2000, step=:exact, log=:off, tol=1e-3)
    # in the original algorithm, x/xk represented distribution of flow
    # in this edited version, it represents the distribution of additional flow from our OD trips on top of inferred existing flows

    setup_time = time()

    if log==:on
        println("-------------------------------------")
        println("Network Name: $(ta_data.network_name)")
        println("Method: $method")
        println("Line Search Step: $step")
        println("Maximum Interation Number: $max_iter_no")
        println("Tolerance for AEC: $tol")
        println("Number of processors: ", nprocs())
    end

    # unpacking data from ta_data
    network_name = ta_data.network_name

    number_of_zones = ta_data.number_of_zones
    number_of_nodes = ta_data.number_of_nodes
#     first_thru_node = ta_data.first_thru_node
    number_of_links = ta_data.number_of_links

    start_node = ta_data.start_node
    end_node = ta_data.end_node
    capacity = ta_data.capacity
    link_length = ta_data.link_length

    free_flow_time = ta_data.free_flow_time
    B = ta_data.B
    power = ta_data.power
    speed_limit = ta_data.speed_limit
    toll = ta_data.toll
    link_type = ta_data.link_type
    number_of_zones = ta_data.number_of_zones
    total_od_flow = ta_data.total_od_flow
    travel_demand = ta_data.travel_demand
    od_pairs = ta_data.od_pairs

    toll_factor = ta_data.toll_factor
    distance_factor = ta_data.distance_factor

    best_objective = ta_data.best_objective

    travel_time = ta_data.travel_time

    # preparing a graph
    graph = create_graph(start_node, end_node, number_of_nodes)
    link_dic = sparse(start_node, end_node, collect(1:number_of_links))
#     link_dic = collect(zip(start_node, end_node))
#     link_dic = Dict{Tuple{Int64,Int64}, Int64}(links[i] => i for i=1:number_of_links)

    # calculate flow from seeded times
    if travel_time != zeros(number_of_links)
        fixed_flow = similar(travel_time)
        for i=1:length(fixed_flow)
            if (free_flow_time[i] > travel_time[i]) || (travel_time[i] == 0) || (free_flow_time[i] == 0)
                fixed_flow[i] = 0
            else
#                 println(travel_time[i], free_flow_time[i])
                fixed_flow[i] = ((travel_time[i]/free_flow_time[i] - 1)/B[i])^(1/power[i]) * capacity[i]
            end
        end
    end
#     print(fixed_flow)
#     return

    setup_time = time() - setup_time

    if log==:on
        println("Setup time = $setup_time seconds")
    end

    function BPR(x)
        # travel_time = free_flow_time .* ( 1.0 + B .* (x./capacity).^power )
        # generalized_cost = travel_time + toll_factor *toll + distance_factor * link_length
        # return generalized_cost

        bpr = similar(x)
        for i=1:length(bpr)
            bpr[i] = free_flow_time[i] * ( 1.0 + B[i] * ((x[i] + fixed_flow[i])/capacity[i])^power[i] )
            bpr[i] += toll_factor * toll[i] + distance_factor * link_length[i]
        end
        return bpr
    end

    function objective(x)
        # value = free_flow_time .* ( x + B.* ( x.^(power+1)) ./ (capacity.^power) ./ (power+1))
        # return sum(value)

        sum = 0.0
        for i=1:length(x)
            sum += free_flow_time[i] * ( (x[i] + fixed_flow[i]) + B[i]* ( (x[i] + fixed_flow[i])^(power[i]+1)) / (capacity[i]^power[i]) / (power[i]+1))
            sum += toll_factor *toll[i] + distance_factor * link_length[i]
        end
        return sum
    end

    function gradient(x)
        return BPR(x)
    end

    function hessian(x)
        no_arc = Base.length(start_node)

        h = zeros(no_arc,no_arc)
        h_diag = hessian_diag(x)

        for i=1:no_arc
            h[i,i] = h_diag[i]
        end

        return h

        #Link travel time = free flow time * ( 1 + B * (flow/capacity)^Power ).
    end

    function hessian_diag(x)
        h_diag = Array{Float64}(size(x))
        for i=1:length(x)
            if power[i] >= 1.0
                h_diag[i] = free_flow_time[i] * B[i] * power[i] * ((x[i] + fixed_flow[i])^(power[i]-1)) / (capacity[i]^power[i])
            else
                h_diag[i] = 0 # Some cases, power is zero.
            end
        end
        # h_diag = free_flow_time .* B .* power .* (x.^(power-1)) ./ (capacity.^power)

        return h_diag
        #Link travel time = free flow time * ( 1 + B * (flow/capacity)^Power ).
    end

    function all_or_nothing_single(travel_time)
        state = LightGraphs.DijkstraState{Float64}
        x = zeros(size(start_node))

        for r=1:size(travel_demand)[1]
            # for each origin node r, find shortest paths to all destination nodes
            state = TA_dijkstra_shortest_paths(graph, travel_time, r, start_node, end_node)

            for s=1:size(travel_demand)[2]
                # for each destination node s, find the shortest-path vector
                # load travel demand
                # x = x + travel_demand[r,s] * get_vector(state, r, s, link_dic)
                add_demand_vector!(x, travel_demand[r,s], state, r, s, link_dic)
            end
        end
        return x
    end

    # parallel computing version
    function all_or_nothing_parallel(travel_time)
        state = LightGraphs.DijkstraState{Float64}
        vv = zeros(size(start_node))
        x = zeros(size(start_node))

        x = x + @parallel (+) for r=1:size(travel_demand)[1]
            # for each origin node r, find shortest paths to all destination nodes
            # if there is any travel demand starting from node r.
            vv = zeros(size(start_node))

            if sum(travel_demand, 2)[r] > 0.0
                state = TA_dijkstra_shortest_paths(graph, travel_time, r, start_node, end_node)

                for s=1:size(travel_demand)[2]
                    # for each destination node s, find the shortest-path vector
                    # v = get_vector(state, r, s, start_node, end_node)

                    if travel_demand[r,s] > 0.0
                        # load travel demand
                        # vv = vv + travel_demand[r,s] * get_vector(state, r, s, link_dic)
                        add_demand_vector!(vv, travel_demand[r,s], state, r, s, link_dic)
                    end
                end

            end

            vv
        end

        return x
    end

    function all_or_nothing(travel_time)
        if nprocs() > 1 # if multiple CPU processes are available
            all_or_nothing_parallel(travel_time)
        else
            all_or_nothing_single(travel_time)
            # when nprocs()==1, using @parallel just adds unnecessary setup time. I guess.
        end
    end

    iteration_time = time()

    # Finding a starting feasible solution
    if travel_time == zeros(number_of_links)
        travel_time = BPR(zeros(number_of_links))
    end
    x0 = all_or_nothing(travel_time)

    # Initializing variables
    xk = x0
    tauk = 0.0
    yk_FW = x0
    sk_CFW = yk_FW
    Hk_diag = Array{Float64,1}

    dk_FW = Array{Float64,1}
    dk_bar = Array{Float64,1}
    dk_CFW = Array{Float64,1}
    dk = Array{Float64,1}

    alphak = 0.0
    Nk = 0.0
    Dk = 0.0

    tauk = 0.0
    is_first_iteration = false
    is_second_iteration = false

    sk_BFW = yk_FW
    sk_BFW_old = yk_FW

    dk_bbar = Array{Float64,1}
    muk = Array{Float64,1}
    nuk = Array{Float64,1}
    beta0 = 0.0
    beta1 = 0.0
    beta2 = 0.0

    best_objective = objective(xk)
#     println(best_objective)

    # function fk(tau)
    #     value = objective(xk+tau*dk)
    #     return value
    # end
    #
    # function lower_bound_k(x, xk)
    #     value = objective(xk) + dot( BPR(xk), ( x - xk) )
    # end

    for k=1:max_iter_no
        # Finding yk
        travel_time = BPR(xk)
        yk_FW = all_or_nothing(travel_time)

        # Basic Frank-Wolfe Direction
        dk_FW = yk_FW - xk
        Hk_diag = hessian_diag(xk) # Hk_diag is a diagonal vector of matrix Hk

        # Finding a feasible direction
        if method == :fw # Original Frank-Wolfe
            dk = dk_FW
        elseif method == :cfw # Conjugate Direction F-W
            if k==1 || tauk > 0.999999 # If tauk=1, then start the process all over again.
                sk_CFW = yk_FW
                dk_CFW = sk_CFW - xk
            else
                dk_bar = sk_CFW - xk # sk_CFW from the previous iteration k-1

                Nk = dot( dk_bar, Hk_diag .* dk_FW )
                Dk = dot( dk_bar, Hk_diag .* (dk_FW - dk_bar) )

                delta = 0.0001 # What value should I use?
                # alphak = 0
                if Dk !=0 && 0 <= Nk/Dk <= 1-delta
                    alphak = Nk/Dk
                elseif Dk !=0 && Nk/Dk > 1-delta
                    alphak = 1-delta
                else
                    alphak = 0
                end

                # Generating new sk_CFW and dk_CFW
                sk_CFW = alphak * sk_CFW + (1-alphak) * yk_FW
                dk_CFW = sk_CFW - xk
            end

            # Feasible Direction to Use for CFW
            dk = dk_CFW
        elseif method == :bfw # Bi-Conjugate Direction F-W

            if tauk > 0.999999
                is_first_iteration = true
                is_second_iteration = true
            end

            if k==1 || is_first_iteration       # First Iteration is like FW
                # println("here")
                sk_BFW_old = yk_FW
                dk_BFW = dk_FW
                is_first_iteration = false
            elseif k==2 || is_second_iteration  # Second Iteration is like CFW
                # println("there")
                dk_bar = sk_BFW_old - xk # sk_BFW_old from the previous iteration 1

                Nk = dot( dk_bar, Hk_diag .* dk_FW )
                Dk = dot( dk_bar, Hk_diag .* (dk_FW - dk_bar) )

                delta = 0.0001 # What value should I use?
                # alphak = 0
                if Dk !=0 && 0 <= Nk/Dk <= 1-delta
                    alphak = Nk/Dk
                elseif Dk !=0 && Nk/Dk > 1-delta
                    alphak = 1-delta
                else
                    alphak = 0
                end

                # Generating new sk_BFW and dk_BFW
                sk_BFW = alphak * sk_BFW_old + (1-alphak) * yk_FW
                dk_BFW = sk_BFW - xk

                is_second_iteration = false
            else
                # println("over there $tauk")
                # sk_BFW, tauk is from iteration k-1
                # sk_BFW_old is from iteration k-2

                dk_bar  = sk_BFW - xk
                dk_bbar = tauk * sk_BFW - xk + (1-tauk) * sk_BFW_old

                muk = - dot( dk_bbar, Hk_diag .* dk_FW ) / dot( dk_bbar, Hk_diag .* (sk_BFW_old - sk_BFW) )
                nuk = - dot( dk_bar, Hk_diag .* dk_FW ) / dot( dk_bar, Hk_diag .* dk_bar) + muk*tauk/(1-tauk)

                muk = max(0, muk)
                nuk = max(0, nuk)

                # println(sk_BFW_old-sk_BFW)

                beta0 = 1 / ( 1 + muk + nuk )
                beta1 = nuk * beta0
                beta2 = muk * beta0

                # dk_BFW = beta0 * dk_FW + beta1 * (sk_BFW - xk) + beta2 * (sk_BFW_old - xk)

                sk_BFW_new = beta0 * yk_FW + beta1 * sk_BFW + beta2 * sk_BFW_old
                dk_BFW = sk_BFW_new - xk

                sk_BFW_old = sk_BFW
                sk_BFW = sk_BFW_new

            end

            # Feasible Direction to Use for BFW
            dk = dk_BFW
        else
            error("The type of Frank-Wolfe method is specified incorrectly. Use :fw, :cfw, or :bfw.")
        end
        # dk is now identified.

        if step==:exact
            # Line Search from xk in the direction dk
            optk = optimize(tau -> objective(xk+tau*dk), 0.0, 1.0, GoldenSection())
            tauk = optk.minimizer
        elseif step==:newton
            # Newton step
            tauk = - dot( gradient(xk), dk ) / dot( dk, Hk_diag.*dk )
            tauk = max(0, min(1, tauk))
        end

        obj = objective(xk)

        # Average Excess Cost
        average_excess_cost = ( dot(xk, travel_time) - dot(yk_FW, travel_time) ) / sum(travel_demand)
        if log==:on
            # println("k=$k,\ttauk=$tauk,\tobjective=$(objective(xk)),\taec=$average_excess_cost")
            @printf("k=%4d, tauk=%15.10f, objective=%15f, aec=%15.10f\n", k, tauk, obj, average_excess_cost)
        end

#         rel_gap = ( objective(xk) - best_objective ) / best_objective
#         rel_gap = ( best_objective - obj )
        rel_gap = abs( objective(xk) - best_objective ) / best_objective

        # Convergence Test
#         if average_excess_cost < tol
        if k > min_iter_no && rel_gap < tol
#         if obj-prior_obj < tol
            if log==:on
                println("Exited before reaching max_iter_no")
            end
            break
        end

        best_objective = min(best_objective, obj)

        # Update x
        new_x = xk + tauk*dk
        xk = new_x

        @assert minimum(xk) >= 0
    end

    if log==:on
        println(best_objective)
    end

    iteration_time = time() - iteration_time

    if log==:on
        println("Iteration time = $iteration_time seconds")
    end

    return xk, travel_time, objective(xk), fixed_flow

end

xk, travel_time, obj, fixed_flow = ta_frank_wolfe(dat; log=:off, method=:cfw, tol=1e-9, max_iter_no=2000, min_iter_no=5) # :bfw results in xk<0 error

df = DataFrame(tail=dat.start_node, head=dat.end_node, travel_time=travel_time, xk=xk, fixed_flow=fixed_flow)
CSV.write(output_filename, df)

max_stars_repo_name: rorysedgwick/trayson

max_issues_repo_name: rorysedgwick/trayson

max_forks_repo_name: rorysedgwick/trayson


text:
module Polyopt

export MomentProb, momentprob, momentprob_chordalembedding, solve_mosek

include("polynomial.jl")
include("cliques.jl")

immutable MomentProb{T<:Number}
    order :: Int
    basis :: Array{Poly{Int},1}
    obj   :: Array{T,1}
    mom   :: Array{Any,1}
    eq    :: Array{Any,1}
end

include("mosek.jl")
include("latex.jl")

# compute all monomials of degree 'deg' or less
function monomials_unsorted{T<:Number}(deg::Int, vars::Array{Poly{T},1})
    if length(vars) == 1 return [ vars[1]^i for i=0:deg ] end

    x, tail = vars[1], vars[2:end]
    monoms = monomials_unsorted(deg, tail)

    for i=1:deg
        append!(monoms, [ x^i * m for m=monomials_unsorted(deg-i, tail) ])
    end
    monoms
end

monomials{T<:Number}(deg::Int, vars::Array{Poly{T},1}) = sort!(monomials_unsorted(deg, vars))

function moment(order::Int, syms::Symbols)
    v = monomials(order, variables(syms))
    v*v'
end

function moment(order::Int, syms::Symbols, I::Array{Int})
    v = monomials(order, variables(syms)[I])
    v*v'
end

function moment{T<:Number}(order::Int, syms::Symbols, p::Poly{T})
    v = monomials(order - (p.deg+1) >> 1, variables(syms))
    p*v*v'
end

function moment{T<:Number}(order::Int, syms::Symbols, p::Poly{T}, I::Array{Int})
    v = monomials(order - (p.deg+1) >> 1, variables(syms)[I])
    p*v*v'
end

function indexmap{T<:Number}(v::Array{Poly{T},1})
    Dict(zip([vi.alpha for vi=v],[1:length(v)]))
end

linear_index(imap::Dict{Array{Int,2},Int}, p::Poly) = Int[ imap[p.alpha[i,:]] for i=1:p.m ]

function vectorize{T<:Number}(p::Poly{T}, l::Int, imap::Dict{Array{Int,2},Int})
    a = zeros(T, l)
    a[linear_index(imap,p) ] = p.c
    a
end

function vectorize{T<:Number}(M::Array{Poly{T},2}, l::Int, imap::Dict{Array{Int,2},Int})
    m, n = size(M)
    m == n || throw(DimensionMismatch("Matrix is not square: dimensions are $(size(M))"))
    subi = Int[]
    subj = Int[]
    val  = T[]
    for j=1:m
        v = M[j,j].c
        push!(subi, (j+(j-1)*m)*ones(Int, length(v))...)
        push!(subj, linear_index(imap,M[j,j])...)
        push!(val, v...)
        for i=j+1:m
            I = linear_index(imap,M[i,j])
            v = M[i,j].c
            push!(subi, (i+(j-1)*m)*ones(Int, length(v))...)
            push!(subj, I...)
            push!(val, v...)
        end
    end
    A = sparse(subi, subj, val, m*m, l)
end

function momentprob{S,T,U,V,W}(order::Int, obj::Poly{S}, pineq::Array{Poly{T},1}, peq::Array{Poly{U},1}, ineq::Array{Poly{V},1}, eq::Array{Poly{W},1})
    v = monomials(2*order, variables(obj.syms))
    l = length(v)
    imap = indexmap(v)

    obj.deg <= 2*order || error("obj has degree higher than 2*order")

    p = vectorize(obj, l, imap)
    mom = Array(Any, length(pineq)+length(ineq)+1)
    mom[1] = vectorize(moment(order, obj.syms), l, imap)
    for k=1:length(pineq)
       pineq[k].deg <= 2*order || error("pineq[$(k)] has degree higher than 2*order")
       mom[k+1] = vectorize(moment(order, pineq[k].syms, pineq[k]), l, imap)
    end

    for k=1:length(ineq)
        ineq[k].deg <= 2*order || error("ineq[$(k)] has degree higher than 2*order")
        mom[k+length(pineq)+1] = sparse(vectorize(ineq[k], l, imap)')
    end

    momeq = Array(Any, length(peq)+length(eq))
    for k=1:length(peq)
        peq[k].deg <= 2*order || error("peq[$(k)] has degree higher than 2*order")
        momeq[k] = vectorize(moment(order, peq[k].syms, peq[k]), l, imap)
    end

    for k=1:length(eq)
        eq[k].deg <= 2*order || error("eq[$(k)] has degree higher than 2*order")
        momeq[k+length(peq)] = sparse(vectorize(eq[k], l, imap)')
    end

    MomentProb(order, v,  p, mom, momeq)
end

momentprob{S}(order::Int, obj::Poly{S}) =
    momentprob(order, obj, Poly{Int}[], Poly{Int}[], Poly{Int}[], Poly{Int}[])

momentprob{S,T}(order::Int, obj::Poly{S}, pineq::Array{Poly{T},1}) =
    momentprob(order, obj, pineq, Poly{Int}[], Poly{Int}[], Poly{Int}[])

momentprob{S,T,U}(order::Int, obj::Poly{S}, pineq::Array{Poly{T},1}, peq::Array{Poly{U},1}) =
    momentprob(order, obj, pineq, peq, Poly{Int}[], Poly{Int}[])

function momentprob_chordal{S,T,U,V,W}(order::Int, cliques::Array{Array{Int,1},1}, obj::Poly{S},
                                       pineq::Array{Poly{T},1}, pineq_index::Array{Int,1},
                                       peq::Array{Poly{U},1}, peq_index::Array{Int,1},
                                       ineq::Array{Poly{V},1}, eq::Array{Poly{W},1})

    v = monomials(2*order, variables(obj.syms))
    l = length(v)
    imap = indexmap(v)

    obj.deg <= 2*order || error("obj has degree higher than 2*order")

    p = vectorize(obj, l, imap)
    mom = Array(Any, length(pineq) + length(cliques) + length(ineq))
    for k=1:length(cliques)
        mom[k] = vectorize(moment(order, obj.syms, cliques[k]), l, imap)
    end

    for k=1:length(pineq)
        pineq[k].deg <= 2*order || error("pineq[$(k)] has degree higher than 2*order")
        mom[length(cliques)+k] = vectorize(moment(order, pineq[k].syms, pineq[k], cliques[pineq_index[k]]), l, imap)
    end

    for k=1:length(ineq)
        ineq[k].deg <= 2*order || error("ineq[$(k)] has degree higher than 2*order")
        mom[length(cliques)+length(pineq)+k] = sparse(vectorize(ineq[k], l, imap)')
    end

    momeq = Array(Any, length(peq)+length(eq))
    for k=1:length(peq)
        peq[k].deg <= 2*order || error("peq[$(k)] has degree higher than 2*order")
        momeq[k] = vectorize(moment(order, peq[k].syms, peq[k], cliques[peq_index[k]]), l, imap)
    end

    for k=1:length(eq)
        eq[k].deg <= 2*order || error("eq[$(k)] has degree higher than 2*order")
        momeq[length(peq)+k] = sparse(vectorize(eq[k], l, imap)')
    end

    MomentProb(order, v,  p, mom, momeq)
end

momentprob_chordal{S,T}(order::Int, cliques::Array{Array{Int,1},1}, obj::Poly{S}, pineq::Array{Poly{T},1}, pineq_index::Array{Int,1}) =
   lasserre_chordal(order, cliques, obj, pineq, pineq_index, Poly{Int}[], Array{Int}[], Poly{Int}[], Poly{Int}[])

function correlative_sparsity{S,T}(obj::Poly{S}, p::Array{Poly{T},1})
    A = eye(Int,obj.n,obj.n)

    for k=1:obj.m
        I = find(obj.alpha[k,:])
        A[I,I] = 1
    end

    for j=1:length(p)
        I = find(sum(p[j].alpha,1))
        A[I,I] = 1
    end

    sparse(A)
end

function momentprob_chordalembedding{S,T,U,V,W}(order::Int, obj::Poly{S},
                                                pineq::Array{Poly{T},1}, peq::Array{Poly{U},1},
                                                ineq::Array{Poly{V},1}, eq::Array{Poly{W},1})

    C = correlative_sparsity(obj, [pineq, ineq, peq, eq])
    cliques = chordal_embedding(C)

    momentprob_chordal(order, cliques, obj,
                       pineq, [ clique_index(cliques, find(sum(p.alpha,1))) for p = pineq ],
                       peq,   [ clique_index(cliques, find(sum(p.alpha,1))) for p = peq   ],
                       ineq, eq)
end

momentprob_chordalembedding{S,T,U}(order::Int, obj::Poly{S}, pineq::Array{Poly{T},1}, peq::Array{Poly{U},1}) =
    momentprob_chordalembedding(order, obj, pineq, peq, Poly{Int}[], Poly{Int}[])

momentprob_chordalembedding{S,T}(order::Int, obj::Poly{S}, pineq::Array{Poly{T},1}) =
    momentprob_chordalembedding(order, obj, pineq, Poly{Int}[], Poly{Int}[], Poly{Int}[])

momentprob_chordalembedding{S}(order::Int, obj::Poly{S}) =
    momentprob_chordalembedding(order, obj, Poly{Int}[], Poly{Int}[], Poly{Int}[], Poly{Int}[])

end


text:
include(joinpath("src", "module.jl"))

config_list = [
    [(5, 5), (5, 6), (6, 5), (6, 6)],
    [(4, 4), (4, 7), (7, 4), (7, 7)],
    [(3, 3), (3, 8), (8, 3), (8, 8)],
    [(2, 2), (2, 9), (9, 2), (9, 9)],
    [(1, 1), (1, 10), (10, 1), (10, 10)]
]

# AxelFraud.run_batch(
#     config_list=config_list,
#     batch_name="expanding_square",
#     grid_dims=(10, 10),
#     steps=3000,
#     replicates=300,
#     when=10,
#     write=true
# )

AxelFraud.run_batch(
    config_list=config_list,
    batch_name="test2",
    grid_dims=(10, 10),
    levels=2,
    steps=300,
    replicates=3,
    when=10,
    write=true
)

AxelFraud.run_random(
    batch_name="random",
    grid_dims=(10, 10),
    stubborn_count=10,
    steps=100,
    replicates=2,
    when=10,
    rnd_seed=0,
    write=true
)


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    combinations = vcat([1], repeat([0], n))
    for i ∈ 1:n for j ∈ i + 1:n + 1 combinations[j] += combinations[j - i] end end
    return combinations[end] - 1
end

compute(5)

compute(100)

@benchmark compute(100)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
using Documenter, ThArrays

makedocs(
  modules=[ThArrays],
  sitename="ThArrays",
  pages = [
      "Home"        => "index.md",
      "Tensor"      => "tensor.md",
      "AD"          => "ad.md",
      "TorchScript" => "torchscript.md",
      "Reference"   => "reference.md",
  ],
  format = Documenter.HTML(prettyurls = haskey(ENV, "GITHUB_EVENT_PATH")))

deploydocs(
    repo = "github.com/TuringLang/ThArrays.jl.git",
    devbranch = "master",
    devurl = "dev",
)


text:
include("../src/Cosmology.jl")
using Cosmology
using Base.Test

function test_approx_eq_rtol(va, vb, rtol, astr, bstr)
    diff = maximum(abs(va - vb))
    tol = rtol*max(maximum(abs(va)), maximum(abs(vb)))
    if diff > tol
        sdiff = string("|", astr, " - ", bstr, "| <= ", tol)
        error("assertion failed: ", sdiff,
              "\n  ", astr, " = ", va,
              "\n  ", bstr, " = ", vb,
              "\n  difference = ", diff, " > ", tol)
    end
end

macro test_approx_eq_rtol(a, b, c)
    :(test_approx_eq_rtol($(esc(a)), $(esc(b)), $(esc(c)), $(string(a)), $(string(b))))
end

# values from http://icosmos.co.uk/

dist_rtol = 1e-6
age_rtol = 2e-4

c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0)
@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1651.9145 dist_rtol
@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3303.829 dist_rtol
@test_approx_eq_rtol comoving_volume_gpc3(c,1) 151.0571 dist_rtol
@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6607.6579 dist_rtol
@test_approx_eq_rtol distmod(c,1) 44.1002 dist_rtol
@test_approx_eq_rtol age_gyr(c,0) 13.4694 age_rtol
@test_approx_eq_rtol age_gyr(c,1) 5.7527 age_rtol
@test_approx_eq_rtol lookback_time_gyr(c,1) 13.4694-5.7527 age_rtol

c = cosmology(h=0.7, OmegaK=0.1, OmegaM=0.3, OmegaR=0)
@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1619.9588 dist_rtol
@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3209.784 dist_rtol
@test_approx_eq_rtol comoving_volume_gpc3(c,1) 140.0856 dist_rtol
@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6479.8352 dist_rtol
@test_approx_eq_rtol distmod(c,1) 44.0578 dist_rtol
@test_approx_eq_rtol age_gyr(c,0) 13.064 age_rtol
@test_approx_eq_rtol age_gyr(c,1) 5.5466 age_rtol
@test_approx_eq_rtol lookback_time_gyr(c,1) 13.064-5.5466 age_rtol

c = cosmology(h=0.7, OmegaK=-0.1, OmegaM=0.3, OmegaR=0)
@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1686.5272 dist_rtol
@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3408.937 dist_rtol
@test_approx_eq_rtol comoving_volume_gpc3(c,1) 163.8479 dist_rtol
@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6746.1088 dist_rtol
@test_approx_eq_rtol distmod(c,1) 44.1453 dist_rtol
@test_approx_eq_rtol age_gyr(c,0) 13.925 age_rtol
@test_approx_eq_rtol age_gyr(c,1) 5.9868 age_rtol
@test_approx_eq_rtol lookback_time_gyr(c,1) 13.925-5.9868 age_rtol

c = cosmology(h=0.7, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1612.0585 dist_rtol
@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3224.1169 dist_rtol
@test_approx_eq_rtol comoving_volume_gpc3(c,1) 140.3851 dist_rtol
@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6448.2338 dist_rtol
@test_approx_eq_rtol distmod(c,1) 44.0472 dist_rtol
@test_approx_eq_rtol age_gyr(c,0) 13.1915 age_rtol
@test_approx_eq_rtol age_gyr(c,1) 5.6464 age_rtol
@test_approx_eq_rtol lookback_time_gyr(c,1) 13.1915-5.6464 age_rtol

c = cosmology(h=0.7, OmegaK=0.1, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1588.0181 dist_rtol
@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3147.6227 dist_rtol
@test_approx_eq_rtol comoving_volume_gpc3(c,1) 132.0466 dist_rtol
@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6352.0723 dist_rtol
@test_approx_eq_rtol distmod(c,1) 44.0146 dist_rtol
@test_approx_eq_rtol age_gyr(c,0) 12.8488 age_rtol
@test_approx_eq_rtol age_gyr(c,1) 5.4659 age_rtol
@test_approx_eq_rtol lookback_time_gyr(c,1) 12.8488-5.4659 age_rtol

c = cosmology(h=0.7, OmegaK=-0.1, OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1637.5993 dist_rtol
@test_approx_eq_rtol comoving_radial_dist_mpc(c,1) 3307.9932 dist_rtol
@test_approx_eq_rtol comoving_volume_gpc3(c,1) 149.8301 dist_rtol
@test_approx_eq_rtol luminosity_dist_mpc(c,1) 6550.3973 dist_rtol
@test_approx_eq_rtol distmod(c,1) 44.0813 dist_rtol
@test_approx_eq_rtol age_gyr(c,0) 13.5702 age_rtol
@test_approx_eq_rtol age_gyr(c,1) 5.8482 age_rtol
@test_approx_eq_rtol lookback_time_gyr(c,1) 13.5702-5.8482 age_rtol

# Test that FlatLCDM works with non-Float64 (BigFloat in this example)
c = cosmology(h=0.7, OmegaM=big(0.3), OmegaR=0)
@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1651.9145 dist_rtol

# Test that FlatWCDM works with non-Float64 (BigFloat in this example)
c = cosmology(h=big(0.7), OmegaM=0.3, OmegaR=0, w0=-0.9, wa=0.1)
@test_approx_eq_rtol angular_diameter_dist_mpc(c,1) 1612.0585 dist_rtol


text:
"""
    create_surfaces_cylinder(cells, setup)

Create surface triangulation of [inner and] outer cylinders [and ECS].

The ground surface is triangulated first, before the walls are "extruded" and
the top surface is copied from the ground surface.
"""
function create_surfaces_cylinder(cells, setup::Setup)
    @unpack radii, centers = cells
    @unpack ncell, rmin, rmax, height, include_in, in_ratio, ecs_shape, ecs_ratio = setup

    # Choose approximate cylinder side length
    nside = 30
    nside_min = 12

    include_ecs = ecs_shape != "no_ecs"
    nboundary = (2 * include_in + 1 + include_ecs) * ncell + include_ecs
    rmean = (rmin + rmax) / 2

    # Create points on boundary circles
    create_nside(r) = max(nside_min, round(Int, nside * r / rmean))
    create_angles(n) = 2π * collect(0:n-1) / n
    create_circle(r, angles, center) = r * [cos.(angles'); sin.(angles')] .+ center
    create_circle_edges(n) = [(1:n-1)' n; (2:n)' 1]

    centers_columns = reshape([centers[:, i] for i = 1:ncell], (1, ncell))

    # Create points for IN compartments
    if include_in
        radii_in = in_ratio * radii
        nside_in = create_nside.(radii_in)
        angles_in = create_angles.(nside_in)
        circles_in = create_circle.(radii_in, angles_in, centers_columns)
        circle_edges_in = create_circle_edges.(nside_in)
        points_in = hcat(circles_in...)
        edges_in = hcat(circle_edges_in...)
        nprevious = 0
        for i = 1:ncell
            ncurrent = nprevious + nside_in[i]
            edges_in[:, nprevious+1:ncurrent] .+= nprevious
            nprevious = ncurrent
        end
    else
        # Create empty arrays
        nside_in = zeros(Int, 1, 0)
        points_in = zeros(2, 0)
        edges_in = zeros(Int, 2, 0)
    end

    # Create points of OUT compartments
    nside_out = create_nside.(radii)
    angles_out = create_angles.(nside_out)
    circles_out = create_circle.(radii, angles_out, centers_columns)
    circle_edges_out = create_circle_edges.(nside_out)
    points_out = hcat(circles_out...)
    edges_out = hcat(circle_edges_out...)
    nprevious = 0
    for i = 1:ncell
        ncurrent = nprevious + nside_out[i]
        edges_out[:, nprevious+1:ncurrent] .+= nprevious
        nprevious = ncurrent
    end

    if include_ecs
        radii_ecs = radii .+ ecs_ratio .* rmean
        nside_ecs = create_nside.(radii_ecs)
        angles_ecs = create_angles.(nside_ecs)
        circles_ecs = create_circle.(radii_ecs, angles_ecs, centers_columns)
        points_ecs = hcat(circles_ecs...)

        if ecs_shape == "box"
            # Determine bounds of domain
            pmin = minimum(points_ecs, dims = 2)
            pmax = maximum(points_ecs, dims = 2)

            # Define four corners of box ECS and their corresponding edges
            points_ecs = [
                pmin[1] pmax[1] pmax[1] pmin[1]
                pmin[2] pmin[2] pmax[2] pmax[2]
            ]
            edges_ecs = [
                1 2 3 4
                2 3 4 1
            ]
        elseif ecs_shape == "convex_hull"
            # Extract points defining convex hull of ECS
            edges_ecs, points_ecs = convexhull(points_ecs)
            npoint_ecs = size(points_ecs, 2)
        elseif ecs_shape == "tight_wrap"
            error("unimplemented")
            # ashape = alphashape(points_ecs)
            # edges_ecs = boundary(hull)
            # boundary_ecs = points_ecs[:, unique(edges_ecs)]
        end
        nedge_ecs = size(points_ecs, 2)
    else
        # Create empty arrays
        points_ecs = zeros(2, 0)
        edges_ecs = zeros(Int, 2, 0)
        nedge_ecs = zeros(Int, 1, 0)
    end

    # Create global collection of points and edges
    # Adjust numbering accordingly
    npoint_in = size(points_in, 2)
    npoint_out = size(points_out, 2)
    npoint_ecs = size(points_ecs, 2)
    points = [points_in points_out points_ecs]
    edges = [edges_in edges_out .+ npoint_in edges_ecs .+ npoint_in .+ npoint_out]
    npoint = size(points, 2)
    nedge = size(edges, 2)

    # Perform Delaynay triangulation of entire domain
    triin = TriangulateIO()
    triin.pointlist = points
    triin.segmentlist = edges
    triin.segmentmarkerlist = 1:nedge
    triout, = triangulate("pQ", triin)
    triangles = eachcol(triout.trianglelist)
    ntriangle = length(triangles)

    boundary_bounds = cumsum([0 nside_in nside_out nedge_ecs], dims = 2)

    ncell_in = ncell * include_in

    find_boundary(node) =
        findfirst(boundary_bounds[1:end-1] .+ 1 .≤ node .≤ boundary_bounds[2:end])

    # Create boundary numbers
    facetmarkers = zeros(Int, ntriangle)
    for (i, t) ∈ enumerate(triangles)
        b = find_boundary.(t)
        b_in = b .≤ ncell_in
        b_out = ncell_in .< b .≤ ncell_in + ncell
        if all(b_in)
            # Inside IN compartment
            facetmarkers[i] = ncell_in + include_ecs * ncell + b[1]
        elseif any(b_in)
            # Triangle touches IN compartment, but is not inside. It is thus the corresponding OUT compartment.
            facetmarkers[i] = ncell_in + include_ecs * ncell + b[.!b_in][1]
        elseif all(b_out) && b[1] == b[2] == b[3]
            # Triangle lies fully within the same OUT comparment
            facetmarkers[i] = ncell_in + include_ecs * ncell + b[1]
        else
            # Triangle lies in the ECS
            facetmarkers[i] = nboundary
        end
    end

    facets = hcat(triangles...)

    # Copy points, edges and markers two the two 3D planes top and bottom
    z = fill(height / 2, 1, npoint)
    points = [points points; -z z]

    function create_sidefacets(le, ri)
        [
            (le:ri-1)' ri (le+1:ri)'.+npoint le.+npoint
            (le:ri-1)'.+npoint ri.+npoint (le:ri-1)'.+npoint ri.+npoint
            (le+1:ri)' le (le+1:ri)' le
        ]
    end

    sidefacets =
        hcat(create_sidefacets.(boundary_bounds[1:end-1] .+ 1, boundary_bounds[2:end])...)

    sidefacetmarkers_in = 1:ncell_in
    sidefacetmarkers_in =
        vcat([repeat([m, m], nside_in[i]) for (i, m) ∈ enumerate(sidefacetmarkers_in)]...)

    sidefacetmarkers_out = (1+!include_ecs)ncell_in+1:(1+!include_ecs)ncell_in+ncell
    sidefacetmarkers_out =
        vcat([repeat([m, m], nside_out[i]) for (i, m) ∈ enumerate(sidefacetmarkers_out)]...)

    sidefacetmarkers_ecs = fill(nboundary, 2npoint_ecs)

    facets = [facets facets .+ npoint sidefacets]
    facetmarkers = [
        facetmarkers
        facetmarkers
        sidefacetmarkers_in
        sidefacetmarkers_out
        sidefacetmarkers_ecs
    ]

    if include_in
        regions_in = [centers; zeros(1, ncell)]
        regions_out = regions_in + [1; 0; 0] * (radii_in + radii) / 2
    else
        regions_in = zeros(3, 0)
        regions_out = [centers; zeros(1, ncell)]
    end
    if include_ecs
        xmin = argmin(points_ecs[1, :])
        region_ecs = [points_ecs[:, xmin]; 0] + [ecs_ratio * rmean / 2; 0; 0]
    else
        region_ecs = zeros(3, 0)
    end
    regions = [regions_in regions_out region_ecs]

    (; points, facets, facetmarkers, regions)
end


text:
using Base.Threads
using TensorOperations
import Base.+
import Base.*
import LinearAlgebra.dot

"""
Addition of two ttvector
"""
function +(x::ttvector{T},y::ttvector{T}) where T<:Number
    @assert(x.ttv_dims == y.ttv_dims, "Dimensions mismatch!")
    d = length(x.ttv_dims)
    ttv_vec = Array{Array{T,3},1}(undef,d)
    rks = x.ttv_rks + y.ttv_rks
    rks[1] = 1
    rks[d+1] = 1
    #initialize ttv_vec
    @threads for k in 1:d
        ttv_vec[k] = zeros(T,x.ttv_dims[k],rks[k],rks[k+1])
    end
   #first core
    ttv_vec[1][:,:,1:x.ttv_rks[2]] = x.ttv_vec[1]
    ttv_vec[1][:,:,(x.ttv_rks[2]+1):rks[2]] = y.ttv_vec[1]
    #2nd to end-1 cores
    @threads for k in 2:(d-1)
        ttv_vec[k][:,1:x.ttv_rks[k],1:x.ttv_rks[k+1]] = x.ttv_vec[k]
        ttv_vec[k][:,(x.ttv_rks[k]+1):rks[k],(x.ttv_rks[k+1]+1):rks[k+1]] = y.ttv_vec[k]
    end
    #last core
    ttv_vec[d][:,1:x.ttv_rks[d],1] = x.ttv_vec[d]
    ttv_vec[d][:,(x.ttv_rks[d]+1):rks[d],1] = y.ttv_vec[d]
    return ttvector{T}(ttv_vec,x.ttv_dims,rks,zeros(d))
end

"""
Addition of two ttoperators
"""
function +(x::ttoperator{T},y::ttoperator{T}) where T<:Number
    @assert(x.tto_dims == y.tto_dims, DimensionMismatch)
    d = length(x.tto_dims)
    tto_vec = Array{Array{T,4},1}(undef,d)
    rks = x.tto_rks + y.tto_rks
    rks[1] = 1
    rks[d+1] = 1
    #initialize tto_vec
    @threads for k in 1:d
        tto_vec[k] = zeros(T,x.tto_dims[k],x.tto_dims[k],rks[k],rks[k+1])
    end
    #first core
    tto_vec[1][:,:,:,1:x.tto_rks[1+1]] = x.tto_vec[1]
    tto_vec[1][:,:,:,(x.tto_rks[2]+1):rks[2]] = y.tto_vec[1]
    #2nd to end-1 cores
    @threads for k in 2:(d-1)
        tto_vec[k][:,:,1:x.tto_rks[k],1:x.tto_rks[k+1]] = x.tto_vec[k]
        tto_vec[k][:,:,(x.tto_rks[k]+1):rks[k],(x.tto_rks[k+1]+1):rks[k+1]] = y.tto_vec[k]
    end
    #last core
    tto_vec[d][:,:,1:x.tto_rks[d],1] = x.tto_vec[d]
    tto_vec[d][:,:,(x.tto_rks[d]+1):rks[d],1] = y.tto_vec[d]
    return ttoperator{T}(tto_vec,x.tto_dims,rks,zeros(d))
end

#matrix vector multiplication in TT format
function *(A::ttoperator{T},v::ttvector{T}) where T<:Number
    @assert(A.tto_dims==v.ttv_dims, DimensionMismatch)
    d = length(A.tto_dims)
    Y = Array{Array{T,3},1}(undef, d)
    A_rks = A.tto_rks #R_0, ..., R_d
    v_rks = v.ttv_rks #r_0, ..., r_d
    @threads for k in 1:d
		M = zeros(A.tto_dims[k], A_rks[k],v_rks[k], A_rks[k+1],v_rks[k+1])
		@tensor M[a,b,c,d,e] = A.tto_vec[k][a,z,b,d]*v.ttv_vec[k][z,c,e]
        Y[k] = reshape(M, A.tto_dims[k], A_rks[k]*v_rks[k], A_rks[k+1]*v_rks[k+1])
    end
    return ttvector{T}(Y,A.tto_dims,A.tto_rks.*v.ttv_rks,zeros(Integer,d))
end

#matrix matrix multiplication in TT format
function *(A::ttoperator{T},B::ttoperator{T}) where T<:Number
    @assert(A.tto_dims==B.tto_dims, DimensionMismatch)
    d = length(A.tto_dims)
    Y = Array{Array{T,4},1}(undef, d)
    A_rks = A.tto_rks #R_0, ..., R_d
    B_rks = B.tto_rks #r_0, ..., r_d
    @threads for k in 1:d
		M = zeros(A.tto_dims[k],A.tto_dims[k],A_rks[k],B_rks[k],A_rks[k+1],B_rks[k+1])
		@tensor M[a,b,c,d,e,f] = A.tto_vec[k][a,z,c,e]*B.tto_vec[k][z,b,d,f]
        Y[k] = reshape(M, A.tto_dims[k], A.tto_dims[k], A_rks[k]*B_rks[k], A_rks[k+1]*B_rks[k+1])
    end
    return ttoperator{T}(Y,A.tto_dims,A.tto_rks.*B.tto_rks,zeros(Integer,d))
end

#dot returns the dot product of two ttvector
function dot(A::ttvector{T},B::ttvector{T}) where T<:Number
    @assert A.ttv_dims==B.ttv_dims "TT dimensions are not compatible"
    d = length(A.ttv_dims)::Int
    A_rks = A.ttv_rks
    B_rks = B.ttv_rks
	out = zeros(T,maximum(A_rks),maximum(B_rks))
    out[1,1] = convert(T,1.0)
    @inbounds for k in 1:d
        M = @view(out[1:A_rks[k+1],1:B_rks[k+1]])
		@tensoropt((α,β), M[a,b] = A.ttv_vec[k][z,α,a]*B.ttv_vec[k][z,β,b]*out[1:A_rks[k],1:B_rks[k]][α,β]) #size R^A_{k} × R^B_{k}
    end
    return out[1,1]::T
end

"""
`dot_par(x_tt,y_tt)' returns the dot product of `x_tt` and `y_tt` in a parallelized algorithm
"""
function dot_par(A::ttvector{T},B::ttvector{T}) where T<:Number
    @assert A.ttv_dims==B.ttv_dims "TT dimensions are not compatible"
    d = length(A.ttv_dims)
    Y = Array{Array{T,2},1}(undef,d)
    A_rks = A.ttv_rks
    B_rks = B.ttv_rks
	C = zeros(T,maximum(A_rks.*B_rks))
    @threads for k in 1:d
		M = zeros(T,A_rks[k],B_rks[k],A_rks[k+1],B_rks[k+1])
		@tensor M[a,b,c,d] = A.ttv_vec[k][z,a,c]*B.ttv_vec[k][z,b,d] #size R^A_{k-1} ×  R^B_{k-1} × R^A_{k} × R^B_{k}
		Y[k] = reshape(M, A_rks[k]*B_rks[k], A_rks[k+1]*B_rks[k+1])
    end
    @inbounds C[1:length(Y[d])] = Y[d][:]
    for k in d-1:-1:1
        @inbounds C[1:size(Y[k],1)] = Y[k]*C[1:size(Y[k],2)]
    end
    return C[1]::T
end

function *(a::S,A::ttvector) where S<:Number
    i = findfirst(isequal(0),A.ttv_ot)
    T = typejoin(typeof(a),eltype(A))
    X = copy(A.ttv_vec)
    X[i] = a*X[i]
    return ttvector{T}(X,A.ttv_dims,A.ttv_rks,A.ttv_ot)
end

function *(a::S,A::ttoperator) where S<:Number
    i = findfirst(isequal(0),A.tto_ot)
    T = typejoin(typeof(a),eltype(A))
    X = copy(A.tto_vec)
    X[i] = a*X[i]
    return ttoperator{T}(X,A.tto_dims,A.tto_rks,A.tto_ot)
end

max_stars_repo_name: msdupuy/Tensor-Train-Julia

max_issues_repo_name: msdupuy/Tensor-Train-Julia

max_forks_repo_name: msdupuy/Tensor-Train-Julia


text:
book = Dict("apple"=> 0.67, "milk"=> 1.49, "avocado"=> 1.49)
println(book)


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(Tuple{typeof(sdf),Matrix{Bool}})   # time: 0.17639625
    @assert precompile(Tuple{typeof(edf_sq),Matrix{Bool}})   # in case the above doesn't work
    @assert precompile(Tuple{typeof(downsample),Matrix{Float64},Int})   # time: 0.0495947
end


text:
###############################################################################
#
#  Isogenies
#
#
#  (C) 2022 Jeroen Hanselman
#
###############################################################################

export Isogeny

export isogeny_from_kernel, isogeny_from_kernel_factored, degree, image,
rational_maps, frobenius_map, isogeny_map_psi, isogeny_map_psi_squared, isogeny_map_phi,
isogeny_map_omega, push_through_isogeny, dual_isogeny, identity_isogeny, multiplication_by_m_map

mutable struct Isogeny{T} <: Map{EllCrv, EllCrv, HeckeMap, Isogeny} where T<: RingElem
  header::MapHeader{EllCrv{T}, EllCrv{T}}
  domain::EllCrv{T}
  codomain::EllCrv{T}
  coordx::RingElem
  coordy::RingElem
  phi::RingElem
  omega::RingElem
  psi::RingElem
  degree::Int

  function Isogeny(E::EllCrv{T}) where T<:RingElem
    r = new{T}()
    r.domain = E
    return r
  end
end

function Isogeny(E::EllCrv{T}, psi::RingElem) where T<:RingElem
  r = Isogeny(E)

  #Normalize and set kernel polynomial
  R = parent(psi)
  psi = numerator(psi//R(leading_coefficient(psi)))
  r.psi = psi

  #Compute codomain and maps
  r.codomain, phi, omega, r.degree = isogeny_kernel_polynomial(E, psi)
  Ry = parent(omega)
  r.coordx = phi//psi^2
  r.coordy = Ry(omega)//Ry(psi^3)

  #Set header
  r.header = MapHeader(E, r.codomain)
  return r
end

#Currently only works for kernels do not contain 4-torsion. Kernels need to be checked for correctness. Input polynomial needs to be separable.
@doc Markdown.doc"""
    isogeny_from_kernel(E::EllCrv, psi::RingElem) -> Isogeny

Compute the isogeny $\phi$: $E$ -> $E'$ where the kernel of $\phi$ contains
the points whose $x$-coordinates are roots of the input polynomial $\psi$.
"""
function isogeny_from_kernel(E::EllCrv, psi::RingElem)
  return Isogeny(E, psi)
end

#Input polynomial needs to be separable. (2^r-torsion is allowed however)
@doc Markdown.doc"""
    isogeny_from_kernel_factored(E::EllCrv, psi::RingElem) -> Isogeny

Compute the isogeny $\phi$: $E$ -> $E'$ where the kernel of $\phi$ contains
the points whose $x$-coordinates are roots of the input polynomial $\psi$.
Return an array of isogenies whose composition equals the isogeny with the given
kernel. The factorisation is determined in the following way: The first maps
will purely contain 2-torsion points in the kernel. When there is no 2-torsion left
the remaining maps will consist of an isogeny with kernel contained in the p-torsion
for every prime divisor p of the degree of the isogeny.
"""
function isogeny_from_kernel_factored(E::EllCrv, psi::RingElem)
  isogeny_list = Isogeny[]

  kernel = psi
  domain = E

  #First compute all isogenies containing 2-torsion
  tor_2 = division_polynomial_univariate(domain, 2)[1]
  ker_G = gcd(kernel, tor_2)

  while (0 != degree(ker_G))

    phi = Isogeny(domain, ker_G)
    push!(isogeny_list, phi)
    kernel = push_through_isogeny(phi, kernel)
    domain = codomain(phi)

    tor_2 = division_polynomial_univariate(domain, 2)[1]
    ker_G = gcd(kernel, tor_2)
  end

  #Then the other isogenies
  d = 2*degree(kernel)+1
  L = sort(prime_divisors(d))

  while (!isempty(L))
    n = L[1]
    tor_n = division_polynomial_univariate(domain, n)[1]
    ker_G = gcd(kernel, tor_n)
    phi = Isogeny(domain, ker_G)
    push!(isogeny_list, phi)
    kernel = push_through_isogeny(phi, kernel)
    domain = codomain(phi)
    d = 2*degree(kernel)+1
    L = sort(prime_divisors(d))
  end
  return isogeny_list
end

@doc Markdown.doc"""
    degree(f::Isogeny) -> Int

Return the degree of the isogeny $f$.
"""
function degree(f::Isogeny)
  return f.degree
end

@doc Markdown.doc"""
    isogeny_map_psi(f::Isogeny) -> Poly

Return the kernel polynomial defining the isogeny $f$.
"""
function isogeny_map_psi(f::Isogeny)
  return f.psi
end

@doc Markdown.doc"""
    isogeny_map_psi_squared(f::Isogeny) -> Poly

Return the denominator of the first coordinate of $f$,
which will usually be the square of the kernel polynomial
(unless the kernel contains 2-torsion, in which case it will almost be a square)
"""
function isogeny_map_psi_squared(f::Isogeny)
  return denominator(f.coordx)
end

@doc Markdown.doc"""
    isogeny_map_phi(f::Isogeny) -> Poly

Return the polynomial phi defining the isogeny $f$.
"""
function isogeny_map_phi(f::Isogeny)
  return numerator(f.coordx)
end

@doc Markdown.doc"""
    isogeny_map_omega(f::Isogeny) -> Poly

Return the polynomial omega defining the isogeny $f$.
"""
function isogeny_map_omega(f::Isogeny)
  return numerator(f.coordy)
end

@doc Markdown.doc"""
    image(phi::Isogeny, P::EllCrvPt) -> EllCrvPt

Return $\phi(P)$.
"""
function image(f::Isogeny, P::EllCrvPt)
  @assert domain(f) == parent(P)
  x = P.coordx
  y = P.coordy

  if evaluate(isogeny_map_psi(f),x) == 0
    return infinity(codomain(f))
  end

  coordx = f.coordx
  coordy = f.coordy

  phix = evaluate(coordx, x)
  phiy = evaluate(evaluate(coordy,y), x)

  return codomain(f)([phix, phiy])
end

@doc Markdown.doc"""
    image(fs::Vector{Isogeny}, P::EllCrvPt) -> EllCrvPt

Return phi_n \circ phi_(n-1) \circ phi_1(P) where fs is a list of compatible
isogenies [phi_1, ..., phi_n].
"""
function image(fs::Vector{Isogeny}, P::EllCrvPt)

  for f in fs
    @assert domain(f) == parent(P)
    P = image(f, P)
  end
  return P
end

@doc Markdown.doc"""
    push_through_isogeny(f::Isogeny, v::RingElem) -> RingElem

Let $f:E \to E'$ be an isogeny and let S be the set of points on E
whose x-coordinate is a root of v. Assume that the kernel polynomial of f
divides v. Return the polynomial psi whose roots are exactly the x-coordinates
of the points $Q$ in $f(S)$.
"""
function push_through_isogeny(f::Isogeny, v::RingElem)

  #The kernel polynomial of f needs to divide v
  R = parent(isogeny_map_psi(f))
  v = divexact(v(gen(R)), isogeny_map_psi(f))

  phi = isogeny_map_phi(f)
  psi_sq = isogeny_map_psi_squared(f)
  Rxy, (x,y) = PolynomialRing(base_ring(phi), 2)
  pol1 = phi(x) - psi_sq(x)*y
  pol2 = v(x)
  L = factor(resultant(pol1, pol2, 1))
  return prod([f for (f,e) in L], init = one(Rxy))(0,gen(parent(phi)))
end

#TODO Need check that we don't need to compose with an automorphism to get the actual dual. Currently we will get the dual up
#to automorphism. Also need to carefully see what happens when the curve is supersingular and we compute the dual of frobenius
@doc Markdown.doc"""
    dual_isogeny(f::Isogeny) -> Isogeny

Return the dual isogeny of f. Currently only returns the dual up to automorphism.
"""
function dual_isogeny(f::Isogeny)

  d = degree(f)
  psi_d = division_polynomial_univariate(f.domain, d)[1]
  psinew = push_through_isogeny(f, psi_d)
  psihat = isogeny_from_kernel(codomain(f), psinew)

  trans = isomorphism(codomain(psihat), domain(f))

  return psihat * trans
end

@doc Markdown.doc"""
    identiy_isogeny((E::EllCrv) -> Isogeny

Return the isogeny corresponding to the identity map on $E$
"""
function identity_isogeny(E::EllCrv)
  return isomorphism_to_isogeny(identity_map(E))
end

@doc Markdown.doc"""
    multiplication_by_m_map((E::EllCrv, m::Int) -> Isogeny

Return the isogeny corresponding to the multiplication by m map on $E$
"""
function multiplication_by_m_map(E::EllCrv, m::S) where S<:Union{Integer, fmpz}

  if m==1
    return isomorphism_to_isogeny(identity_map(E))
  end

  a1, a2, a3 = a_invars(E)

  Kx, x = PolynomialRing(base_field(E), "x")
  Kxy, y = PolynomialRing(Kx, "y")

  mult_mx = multiplication_by_m_numerator(E, m, x)//multiplication_by_m_denominator(E, m, x)
  mult_mx_deriv = derivative(mult_mx)

  x = Kxy(x)
  mult_mxy = multiplication_by_m_numerator(E, m, x)//multiplication_by_m_denominator(E, m, x)
  mult_mx_derivy = numerator(mult_mx_deriv)(x)//denominator(mult_mx_deriv)(x)

  mult_my = ((2*y+a1*x+a3)*mult_mx_derivy//m - a1*mult_mxy-a3)//2

  mul_m = Isogeny(E)
  mul_m.coordx = mult_mx
  mul_m.coordy = mult_my
  mul_m.degree = m^2

  #This feels superfluous. Maybe division polynomials need to be cached somehow
  mul_m.psi = division_polynomial_univariate(E, m)[1]
  mul_m.codomain = E
  mul_m.header = MapHeader(E, E)

  return mul_m
end

@doc Markdown.doc"""
    frobenius_map(E::EllCrv{FinFieldElem}) -> Isogeny

Return the Frobenius map on E.
"""
function frobenius_map(E::EllCrv{T}) where T<:FinFieldElem
  return frobenius_map(E, 1)
end

@doc Markdown.doc"""
    frobenius_map(E::EllCrv{FinFieldElem}, n::Int) -> Isogeny

Return the $n$-th power of the Frobenius map on E.
"""
function frobenius_map(E::EllCrv{T}, n) where T<:FinFieldElem
  f = Isogeny(E)
  K = base_field(E)
  p = characteristic(K)
  pn = p^n
  Rx, x = PolynomialRing(K, "x")
  Rxy, y = PolynomialRing(Rx, "y")
  f.codomain = E
  f.degree = pn
  f.coordx = x^pn//1
  f.coordy = y^pn//1
  f.psi = one(Rx)
  f.header = MapHeader(E, f.codomain)
  return f
end

@doc Markdown.doc"""
    rational_maps(I::Isogeny) -> [Poly, Poly]

Return the rational maps defining the isogeny.
"""
function rational_maps(I::Isogeny)
  return [I.coordx, I.coordy]
end

@doc Markdown.doc"""
    compose(I1::Isogeny, I2::Isogeny) -> Isogeny

Return composition $I2 \circ I1$.
"""
function compose(I1::Isogeny, I2::Isogeny)
#Maybe this can be done better. Currently I just evaluate and recalculate phi, omega and psi based on numerator and
#denominator

  @assert codomain(I1) == domain(I2)

  E = domain(I1)
  F = codomain(I2)

  newx = I1.coordx
  newy = I1.coordy
  Rxy = parent(newy)

  tempx = evaluate(I2.coordx, newx)

  newx_overy = Rxy(numerator(newx))//Rxy(denominator(newx))
  tempomega_num = numerator(I2.coordy)
  tempomega_denom = denominator(I2.coordy)

  omega_num = evaluate(tempomega_num.coeffs[1], newx_overy) + evaluate(tempomega_num.coeffs[2], newx_overy)*newy
  omega_denom = evaluate(tempomega_denom.coeffs[1], newx_overy)

  #To compute the kernel polynomial psi, we need to factor out the 2-torsion part first before we can take the square of the denominator

  psi_2 = division_polynomial_univariate(E, 2)[1]

  R = parent(psi_2)
  psitemp = denominator(tempx)(gen(R))
  psi_2 = numerator(psi_2//R(leading_coefficient(psi_2)))
  psi_G = gcd(psi_2, psitemp)
  if (0 != degree(psi_G))
    newpsi = psi_G*sqrt(divexact(psitemp, psi_G))
  else
    newpsi = sqrt(psitemp)
  end

  E = domain(I1)
  F = codomain(I2)

  f = Isogeny(E)
  f.coordx = evaluate(I2.coordx, newx)
  f.coordy = omega_num//omega_denom
  f.psi = newpsi

  f.codomain = F
  f.degree = I1.degree*I2.degree
  f.header = MapHeader(E, F)
  return f
end

@doc Markdown.doc"""
    compose(fs::Vector{Isogeny}) -> Isogeny

Return the composition phi_n \circ phi_(n-1) \circ phi_1 where fs is a list of compatible
isogenies [phi_1, ..., phi_n].
"""
function compose(fs::Vector{Isogeny})
  g = fs[1]
  n = length(fs)
  for i in (2:n)
    g = g*fs[i]
  end
  return g
end

@doc Markdown.doc"""
    is_isogenous(E::EllCrv, F::EllCrv) -> Bool

Return true when $E$ and $F$ are isogenous. Currently only implemented for
curves over finite fields.
"""
function is_isogenous(E::EllCrv{T}, F::EllCrv{T}) where T<:FinFieldElem
  return order(E) == order(F)
end

#Algorithm uses Kohel's formulas for isogenies. Based on implementation in Sage described
#in Daniel Shumow's thesis: Isogenies of Elliptic Curves: A Computational Approach
function isogeny_kernel_polynomial(E::EllCrv, psi)

#TODO: Test if polynomial defines a possible kernel, i.e. it needs to be monic, separable, defined over the base field and its roots (in the alg. closure)
#need to be x-coordinates of torsion points on the curve

  R = parent(psi)

  psi_G = gcd(division_polynomial_univariate(E, 2, gen(R))[1], psi)
  # force this polynomial to be monic
  psi_G = numerator(psi_G//R(leading_coefficient(psi_G)))

  if (0 != degree(psi_G))
    psi_quo = numerator((psi//psi_G))
      if (degree(psi_quo)!=0)
          throw(DomainError(psi, "In order to use Kohel's algorithm, a kernel of even degree needs to be contained in the two torsion."))
      end
    phi, omega, v, w, n, d = even_kernel_polynomial(E,psi)
  else
    phi, omega, v, w, n, d = odd_kernel_polynomial(E,psi)
  end
  phi = phi
  psi = psi
  return compute_codomain(E,v,w), phi, omega, d
end

function even_kernel_polynomial(E::EllCrv, psi_G)
  R = parent(psi_G)
  x = gen(R)
  psi_2tor = division_polynomial_univariate(E,2, x)[1]
  if  mod(psi_2tor,psi_G) != 0
    throw(DomainError(psi_G,"Does not define a finite subgroup of the elliptic curve."))
  end
  n = degree(psi_G)
  d = n+1
  K = base_field(E)
  char = characteristic(K)

  Kxy,y = PolynomialRing(R,"y")

  a1, a2, a3, a4, a6 = a_invars(E)
  b2, b4, b6 = E.b_invars

  if n == 1
    x0 = constant_coefficient(-psi_G)

    # determine y0
    if (2 == characteristic(base_field(E)))
      y0 = sqrt(x0^3 + a2*x0^2 + a4*x0 + a6)
    else
      y0 = -(a1*x0 + a3)//2
      v = (3*x0^2 + 2*a2*x0 + a4 - a1*y0)
      w = x0*v

      phi = (x*psi_G + v)*psi_G
      omega = (y*psi_G^2 - v*(a1*psi_G + (y - y0)))*psi_G
    end
  elseif (3 == n)
    s = coefficients(psi_G)
    s1 = -s[n-1]
    s2 = s[n-2]
    s3 = -s[n-3]

    dpsi_G = derivative(psi_G)
    ddpsi_G = derivative(dpsi_G)

    phi = (dpsi_G^2) + (-2*ddpsi_G + (4*x - s1))*psi_G
    dphi = derivative(phi)

    psi_2 = 2*y + a1*x + a3

    omega = (psi_2*(dphi*psi_G - phi*dpsi_G) - (a1*phi + a3*psi_G)*psi_G)//2

    phi = phi*psi_G
    omega = omega*psi_G

    v = 3*(s1^2 - 2*s2) + b2*s1//2 + 3*b4//2
    w = 3*(s1^3 - 3*s1*s2 + 3*s3) + b2*(s1^2 - 2*s2)//2 + b4*s1//2
  else
     throw(DomainError(psi_G, "Input polynomial must be of degree 1 or 3."))
  end
return phi, omega, v, w, n, d
end

function odd_kernel_polynomial(E, psi)
  n = degree(psi)
  d = 2*n+1

  char = characteristic(base_field(E))

  a1, a2, a3, a4, a6 = a_invars(E)
  b2, b4, b6 = b_invars(E)

  R = parent(psi)
  x = gen(R)

  Rxy,y = PolynomialRing(R,"y")
  psi_2 = 2*y + a1*x + a3

  psicoeffs = coefficients(psi)

  s1 = s2 = s3 = 0

  if (1 <= n)
    s1 = -psicoeffs[n-1]
  end
  if (2 <= n)
    s2 = psicoeffs[n-2]
  end
  if (3 <= n)
    s3 = -psicoeffs[n-3]
  end

  #Variables for computing the codomain
  v = 6*(s1^2 - 2*s2) + b2*s1 + n*b4
  w = 10*(s1^3 - 3*s1*s2 + 3*s3) + 2*b2*(s1^2 - 2*s2) + 3*b4*s1 + n*b6

  dpsi = derivative(psi)
  ddpsi = derivative(dpsi)

  phi = (4*x^3 + b2*x^2 + 2*b4*x + b6)*(dpsi^2 - ddpsi*psi) - (6*x^2 + b2*x + b4)*dpsi*psi + (d*x - 2*s1)*psi^2
  dphi = derivative(phi)

  if (char!=2)
    omega = numerator(dphi*psi*psi_2//2) - phi*dpsi*psi_2 -numerator((a1*phi + a3*psi^2)*psi//2)
  else
  #Characteristic 2 case:
    if (0 < n)
      s1 = -psicoeffs[n-1]
    else
      s1 = 0
    end
    ddpsi = 0
    temp = 1
    for j in (0:n-1)
      ddpsi = ddpsi + binomial(j+2,2)*psicoeffs[j+2]*temp
      temp = x*temp
    end

    dddpsi = 0
    temp = 1
    for j in (0:n-2)
      dddpsi = dddpsi + (3*binomial(j+3,3))*psicoeffs[(j+3)]*temp
      temp = x*temp
    end

    omega = dphi*psi*y - phi*dpsi*psi_2 + ((a1*x + a3)*(psi_2^2)*(ddpsi*dpsi-dddpsi*psi) + (a1*psi_2^2 - 3*(a1*x + a3)*(6*x^2 + b2*x + b4))*dpsi*psi +(a1*x^3 + 3*a3*x^2 + (2*a2*a3 - a1*a4)*x + (a3*a4 - 2*a1*a6))*dpsi^2 +(-(3*a1*x^2 + 6*a3*x + (-a1*a4 + 2*a2*a3)) + (a1*x + a3)*(d*x - 2*s1) )*dpsi*psi + (a1*s1 + a3*n)*psi^2)*psi
  end
  return phi, omega, v, w, n, d
end

function compute_codomain(E::EllCrv, v, w)
  a1, a2, a3, a4, a6 = a_invars(E)
  newa4 = a4 - 5*v
  newa6 = a6 - (a1^2 + 4*a2)*v - 7*w
  return EllipticCurve([a1, a2, a3, newa4, newa6])
end