func (mp MatrixProfile) distanceProfile(idx int, profile []float64, fft *fourier.FFT) error {
	if idx > len(mp.A)-mp.W {
		return fmt.Errorf("provided index  %d is beyond the length of timeseries %d minus the subsequence length %d", idx, len(mp.A), mp.W)
	}
	if err := mp.mass(mp.A[idx:idx+mp.W], profile, fft); err != nil {
		return err
	}
	if mp.SelfJoin {
		util.ApplyExclusionZone(profile, idx, mp.W/2)
	}
	return nil
}
distanceProfile computes the distance profile between a and b time series.
If b is set to nil then it assumes a self join and will create an exclusion
area for trivial nearest neighbors. Writes the euclidean distance between
the specified subsequence in mp.A with each subsequence in mp.B to profile
package matrixprofile
import (
	"container/heap"
	"encoding/json"
	"errors"
	"fmt"
	"io/ioutil"
	"math"
	"math/rand"
	"os"
	"runtime"
	"sort"
	"strconv"
	"sync"
	"github.com/matrix-profile-foundation/go-matrixprofile/av"
	"github.com/matrix-profile-foundation/go-matrixprofile/util"
	"gonum.org/v1/gonum/dsp/fourier"
	"gonum.org/v1/gonum/floats"
	"gonum.org/v1/plot/plotter"
)
type MatrixProfile struct {
	A        []float64    `json:"a"`                  
	B        []float64    `json:"b"`                  
	AMean    []float64    `json:"a_mean"`             
	AStd     []float64    `json:"a_std"`              
	BMean    []float64    `json:"b_mean"`             
	BStd     []float64    `json:"b_std"`              
	BF       []complex128 `json:"b_fft"`              
	N        int          `json:"n"`                  
	W        int          `json:"w"`                  
	SelfJoin bool         `json:"self_join"`          
	MP       []float64    `json:"mp"`                 
	Idx      []int        `json:"pi"`                 
	MPB      []float64    `json:"mp_ba"`              
	IdxB     []int        `json:"pi_ba"`              
	AV       av.AV        `json:"annotation_vector"`  
	Opts     *MPOpts      `json:"options"`            
	Motifs   []MotifGroup
	Discords []int
}
func New(a, b []float64, w int) (*MatrixProfile, error) {
	if a == nil || len(a) == 0 {
		return nil, fmt.Errorf("first slice is nil or has a length of 0")
	}
	if b != nil && len(b) == 0 {
		return nil, fmt.Errorf("second slice must be nil for self-join operation or have a length greater than 0")
	}
	mp := MatrixProfile{
		A: a,
		W: w,
		N: len(b),
	}
	if b == nil {
		mp.N = len(a)
		mp.B = a
		mp.SelfJoin = true
	} else {
		mp.B = b
	}
	if mp.W > len(mp.A) || mp.W > len(mp.B) {
		return nil, fmt.Errorf("subsequence length must be less than the timeseries")
	}
	if mp.W < 2 {
		return nil, fmt.Errorf("subsequence length must be at least 2")
	}
	mp.AV = av.Default
	return &mp, nil
}
func applySingleAV(mp, ts []float64, w int, a av.AV) ([]float64, error) {
	avec, err := av.Create(a, ts, w)
	if err != nil {
		return nil, err
	}
	if len(avec) != len(mp) {
		return nil, fmt.Errorf("annotation vector length, %d, does not match matrix profile length, %d", len(avec), len(mp))
	}
	maxMP := 0.0
	for _, val := range mp {
		if val > maxMP {
			maxMP = val
		}
	}
	for idx, val := range avec {
		if val < 0.0 || val > 1.0 {
			return nil, fmt.Errorf("got an annotation vector value of %.3f at index %d. must be between 0 and 1", val, idx)
		}
	}
	out := make([]float64, len(mp))
	for idx, val := range avec {
		out[idx] = mp[idx] + (1-val)*maxMP
	}
	return out, nil
}
func (mp MatrixProfile) ApplyAV() ([]float64, []float64, error) {
	var err error
	abmp := make([]float64, len(mp.MP))
	bamp := make([]float64, len(mp.MPB))
	copy(abmp, mp.MP)
	copy(bamp, mp.MPB)
	if !mp.Opts.Euclidean {
		util.P2E(abmp, mp.W)
		util.P2E(bamp, mp.W)
	}
	abmp, err = applySingleAV(abmp, mp.A, mp.W, mp.AV)
	if err != nil {
		return nil, nil, err
	}
	if mp.MPB != nil {
		bamp, err = applySingleAV(bamp, mp.B, mp.W, mp.AV)
	}
	if err != nil {
		return nil, nil, err
	}
	if !mp.Opts.Euclidean {
		util.E2P(abmp, mp.W)
		util.E2P(bamp, mp.W)
	}
	return abmp, bamp, nil
}
func (mp MatrixProfile) Save(filepath, format string) error {
	var err error
	switch format {
	case "json":
		f, err := os.Open(filepath)
		if err != nil {
			f, err = os.Create(filepath)
			if err != nil {
				return err
			}
		}
		defer f.Close()
		out, err := json.Marshal(mp)
		if err != nil {
			return err
		}
		_, err = f.Write(out)
	default:
		return fmt.Errorf("invalid save format, %s", format)
	}
	return err
}
func (mp *MatrixProfile) Load(filepath, format string) error {
	var err error
	switch format {
	case "json":
		f, err := os.Open(filepath)
		if err != nil {
			return err
		}
		defer f.Close()
		b, err := ioutil.ReadAll(f)
		if err != nil {
			return err
		}
		err = json.Unmarshal(b, mp)
	default:
		return fmt.Errorf("invalid load format, %s", format)
	}
	return err
}
type mpVals []float64
func (m mpVals) Len() int {
	return len(m)
}
func (m mpVals) Swap(i, j int) {
	m[i], m[j] = m[j], m[i]
}
func (m mpVals) Less(i, j int) bool {
	return m[i] < m[j]
}
func (m *mpVals) Push(x interface{}) {
	*m = append(*m, x.(float64))
}
func (m *mpVals) Pop() interface{} {
	x := (*m)[len(*m)-1]
	*m = (*m)[:len(*m)-1]
	return x
}
type MPDistOpts struct {
	AV   av.AV
	Opts *MPOpts
}
func NewMPDistOpts() *MPDistOpts {
	return &MPDistOpts{
		AV:   av.Default,
		Opts: NewMPOpts(),
	}
}
func MPDist(a, b []float64, w int, o *MPDistOpts) (float64, error) {
	if o == nil {
		o = NewMPDistOpts()
	}
	mp, err := New(a, b, w)
	if err != nil {
		return 0, err
	}
	if err = mp.Compute(o.Opts); err != nil {
		return 0, nil
	}
	mpab, mpba, err := mp.ApplyAV()
	if err != nil {
		return 0, nil
	}
	thresh := 0.05
	k := int(thresh * float64(len(a)+len(b)))
	mpABBASize := len(mpab) + len(mpba)
	if k < mpABBASize {
		var lowestMPs mpVals
		heap.Init(&lowestMPs)
		for _, d := range mpab {
			if !mp.Opts.Euclidean {
				d = -d
			}
			if len(lowestMPs) == k+1 {
				if d < lowestMPs[0] {
					heap.Pop(&lowestMPs)
					heap.Push(&lowestMPs, d)
				}
			} else {
				heap.Push(&lowestMPs, d)
			}
		}
		for _, d := range mpba {
			if !mp.Opts.Euclidean {
				d = -d
			}
			if len(lowestMPs) == k+1 {
				if d < lowestMPs[0] {
					heap.Pop(&lowestMPs)
					heap.Push(&lowestMPs, d)
				}
			} else {
				heap.Push(&lowestMPs, d)
			}
		}
		if !mp.Opts.Euclidean {
			return -lowestMPs[0], nil
		}
		return lowestMPs[0], nil
	}
	var trackVal float64
	if !mp.Opts.Euclidean {
		trackVal = 1
	}
	for _, d := range mp.MP {
		if mp.Opts.Euclidean {
			if d > trackVal {
				trackVal = d
			}
		} else {
			if d < trackVal {
				trackVal = d
			}
		}
	}
	for _, d := range mp.MPB {
		if mp.Opts.Euclidean {
			if d > trackVal {
				trackVal = d
			}
		} else {
			if d < trackVal {
				trackVal = d
			}
		}
	}
	return trackVal, nil
}
type Algo string
const (
	AlgoSTOMP Algo = "stomp"
	AlgoSTAMP Algo = "stamp"
	AlgoSTMP  Algo = "stmp"
	AlgoMPX   Algo = "mpx"
)
type MPOpts struct {
	Algorithm    Algo    `json:"algorithm"`   
	SamplePct    float64 `json:"sample_pct"`  
	NJobs        int     `json:"n_jobs"`
	Euclidean    bool    `json:"euclidean"`                   
	RemapNegCorr bool    `json:"remap_negative_correlation"`  
}
func NewMPOpts() *MPOpts {
	p := runtime.NumCPU() * 2
	if p < 1 {
		p = 1
	}
	return &MPOpts{
		Algorithm: AlgoMPX,
		SamplePct: 1.0,
		NJobs:     p,
		Euclidean: true,
	}
}
func (mp *MatrixProfile) Compute(o *MPOpts) error {
	if o == nil {
		o = NewMPOpts()
	}
	mp.Opts = o
	if o.SamplePct < 1 {
		return mp.stamp()
	}
	switch o.Algorithm {
	case AlgoSTOMP:
		return mp.stomp()
	case AlgoSTAMP:
		return mp.stamp()
	case AlgoSTMP:
		return mp.stmp()
	case AlgoMPX:
		return mp.mpx()
	default:
		return fmt.Errorf("Unsupported algorithm for matrix profile, %s", o.Algorithm)
	}
	return nil
}
func (mp *MatrixProfile) initCaches() error {
	var err error
	mp.BMean, mp.BStd, err = util.MovMeanStd(mp.B, mp.W)
	if err != nil {
		return err
	}
	mp.AMean, mp.AStd, err = util.MovMeanStd(mp.A, mp.W)
	if err != nil {
		return err
	}
	fft := fourier.NewFFT(mp.N)
	mp.BF = fft.Coefficients(nil, mp.B)
	return nil
}
func (mp MatrixProfile) crossCorrelate(q []float64, fft *fourier.FFT) []float64 {
	qpad := make([]float64, mp.N)
	for i := 0; i < len(q); i++ {
		qpad[i] = q[mp.W-i-1]
	}
	qf := fft.Coefficients(nil, qpad)
	for i := 0; i < len(qf); i++ {
		qf[i] = mp.BF[i] * qf[i]
	}
	dot := fft.Sequence(nil, qf)
	for i := 0; i < mp.N-mp.W+1; i++ {
		dot[mp.W-1+i] = dot[mp.W-1+i] / float64(mp.N)
	}
	return dot[mp.W-1:]
}
func (mp MatrixProfile) mass(q []float64, profile []float64, fft *fourier.FFT) error {
	qnorm, err := util.ZNormalize(q)
	if err != nil {
		return err
	}
	dot := mp.crossCorrelate(qnorm, fft)
	for i := 0; i < len(dot); i++ {
		profile[i] = math.Sqrt(math.Abs(2 * (float64(mp.W) - (dot[i] / mp.BStd[i]))))
	}
	return nil
}