Physics Candles
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © Techotomy

//@version=5
indicator(title="Physics Candles", shorttitle = "Physics Beta", overlay=true, max_labels_count=500, max_boxes_count=500)

EMPTY_COLOR = color.new(color.black, 100)
BUG_PREFIX = "ᄽ(◉෴◉)ᄿ Evil Bug Detected: "
ERROR_PREFIX = "(◎_◎) Error Detected: "

MIN_COLOR_UP = color.new(color.green, 90)
MAX_COLOR_UP = color.new(color.green, 10)
MIN_COLOR_DN = color.new(color.red, 90)
MAX_COLOR_DN = color.new(color.red, 10)

price_source_tooltip = "[Optional] Specify the candle price source by entering the 
 data feed source and ticker symbol bisected by a colon (E.g., NASDAQ:TSLA or 
 CME_MINI:ES!.).  The price source defaults to the TV chart symbol."
physics_source_tooltip = "[Optional] Specify the candle volume source by entering 
 the data feed source and ticker symbol bisected by a colon (E.g., NASDAQ:TSLA or 
 CME_MINI:ES!).  The price source defaults to the TV chart symbol. "
cycle_tooltip = "Specify the cycle time frame that restarts window collection 
 when compares current volume intensity against.  For example, to compare real-time 
 relative volume in the current 15-minute window for each day, select Daily"
sample_tooltip = "Specify the raw sampling granularity.  A lower timeframe may
 give higher accuracy in volume segmentation at the cost of sample count due to 
 TV data limits.  Raw volume samples are processed by segmentation and are not
 equivalent to the processed sample count."
window_tooltip = "Specify the window time frame to compare current relative volume 
 against.  For example, to compare real-time relative volume in the f
 current 15-minute window for each day, select 15-minute.  Low window timeframes 
 are useful for market open and closes to account for their inherent variability
 minute-to-minute"
sample_date_tooltip = "Specify the start and end date for data collection. These
 default to the maximum start and end date times enforced by TV data limits."
distribution_tooltip = "Describe the distribution of volume in candle segments.  
 A normal distribution assumes samples are distributed as a standard bell-curve 
 while a log-normal distribution takes the logarithm of samples to convert to 
 a standard bell-curve."
method_tooltip = "Describe the method of volume data transformation before 
 normalization to the color range.  Linear uses raw volume values without 
 transformation, logarithm takes the log of volume, z-score converts volume to 
 the multiples of the standard deviations away from the mean, and power-law takes 
 the exponent of volume to the power of gamma.  Logarithm and power-law with a 
 gamma below one tend to increase sensitivity while power-law with a gamma above 
 one decreases sensivity to highlight large volume anomalies."
gamma_tooltip = "Specify the value of the exponent when using a power-law
 indicator transformation"
cutoff_tooltip = "Specify the static minimum and maximum values to normalize between the 
 minimum and maximum color range.  The scores at or below the minimum 
 cutoff will have the same intensity as the minimum color settings while the 
 scores at or above the maximum cutoff will have the same intensity as the 
 maximum color settings with a linear transition to color for scores in between 
 the two cutoffs"
autonorm_tooltip = "Auto-Normalization automatically interpolates color intensity 
 between the minimum and maximum multiples of standard deviations from the mean for 
 each distribution/window of physics values.  For example, a minimum multiplier of 0 
 and a maximum multiplier of 3 interpolates colors between the mean (mean + 0*stdev) 
 and 3 standard deviations from the mean (mean + 3*stdev) for each window.  Values 
 below the mean are pinned to the minimum color scheme while values above 3 are pinned 
 to the maximum color scheme. This is different than the static min and max cutoff 
 values above since those values remain constant across every window and are not 
 automatically adjusted according to the aggregrated statistics of each window."
plot_dist_tooltip = "Plot the processed sample distribution for a given window, 
 where samples are segmented and scored.  To find the window of interest to plot,
 use the Print Windows feature."
resolution_tooltip = "Specify the number of volume segments per candle.  The
 maximum number is 8 when using one indicator, but additional indicators can be
 applied to the chart to increase the resolution.  For example, for 16 segments
 per candle, set the resolution to 16 with one indicator Stack Number set to 0
 and the other indicator Stack Number set to 1."
stack_tooltip = "Specify the indicator number when increasing candle resolution
 in increasing order starting with 0 (e.g., 0 for first indicator, 1 for second, 
 2 for the third)"
overlay_tooltip = "Specify the number of indicator overlays to
 automaticaly decrease the transparency of the color scheme for each individual 
 indicator such that combined overlay equals the target transparency.  The maximum
 transparency of the color scheme should be less than 100% since only an infinite 
 number of transparent overlays can equal 100%."
hide_tooltip = "Specify the intraday time to hide  theindicator.  This is useful 
 when calibrating multiple indicators to different parts of the trading day."


price_source = input.string(defval='', title='Candle Source  ',
 group='Sources', tooltip=price_source_tooltip)
physics_source = input.string(defval='', title='Physics Source  ', 
 group='Sources', tooltip=physics_source_tooltip)
mass_selection = input.string(defval='Volume', title='Mass', group='Sources',
 options=['Volume', 'Time'])
physics_selection = input.string(defval='Mass', title='Physics', group='Sources', 
 options= ['Mass', 'Position', 'Velocity', 'Acceleration', 'Momentum', 
 'Force', 'Kinetic Energy'])

input_to_timeframe(str) =>
    string tf = switch str
        "1 Min" => "1"
        "2 Mins" => "2"
        "3 Mins" => "3"
        "5 Mins" => "5"
        "10 Mins" => "10"
        "15 Mins" => "15"
        "30 Mins" => "30"
        "45 Mins" => "45"
        "1 Hour" => "60"
        "2 Hours" => "120"
        "3 Hours" => "180"
        "4 Hours" => "240"
        "24 Hours" => "1440"
        "1 day" => "1D"
        "1 week" => "1W"
        "1 month" => "1M"
        "1 quarter" => "3M"
        "1 year" => "12M"
        "Daily" => "1D"
        "Weekly" => "1W"
        "Monthly" => "1M"
        "Quarterly" => "3M"
        "Yearly" => "12M"
        "Day" => "1D"
        "Week" => "1W"
        "Month" => "1M"
        "Quarter" => "3M"
        "Year" => "12M"

ctf = input.timeframe(
 defval="Daily", title='Cycle Time Frame', group='Data Collection', options=
 ["Daily", "Weekly", "Monthly", "Quarterly", "Yearly"], tooltip=cycle_tooltip)
cycle_tf = input_to_timeframe(ctf)

stf = input.string(
 defval="1 Min", title='Sample Time Frame', group='Data Collection', options=
 ["1 Min", "2 Mins", "3 Mins", "5 Mins", "10 Mins", "15 Mins", "30 Mins", 
 "45 Mins", "1 Hour", "2 Hours", "3 Hours", "4 Hours", "24 Hours", "Daily", 
 "Weekly", "Monthly", "Quarterly", "Yearly"], tooltip=sample_tooltip)
sample_tf = input_to_timeframe(stf)

wtf = input.string(
 defval="30 Mins", title='Window Time Frame', group='Data Collection', options=
 ["1 Min", "2 Mins", "3 Mins", "5 Mins", "10 Mins", "15 Mins", "30 Mins", 
 "45 Mins", "1 Hour", "2 Hours", "3 Hours", "4 Hours", "24 Hours", "Daily", 
 "Weekly", "Monthly", "Quarterly", "Yearly"], tooltip=window_tooltip)
window_tf = input_to_timeframe(wtf)

use_start_date = input.bool(defval=false, title='Start', inline='start', 
 group='Data Collection')
start_date = input.time(
 defval=timestamp("1 Jan 2021 0:00"), title='', inline='start', 
 group='Data Collection', tooltip=sample_date_tooltip)
use_end_date = input.bool(defval=false, title='End ', inline='end', 
 group='Data Collection')
end_date = input.time(
 defval=timestamp("1 Jan 2022 0:00"), title='', inline='end', 
 group='Data Collection', tooltip=sample_date_tooltip)
 
dist_type = input.string(
 defval='Log-Normal', title='Distribution', group='Physics Normalization',
 options=['Normal', 'Log-Normal'], tooltip=distribution_tooltip)
norm_type = input.string(
 defval="Z-Score", title='Method', group='Physics Normalization', 
 options=["Linear", "Logarithm", "Z-Score", "Power-Law"],
 tooltip=method_tooltip)
min_score = input.float(
 defval=0, title='Min Score Cutoff', group='Physics Normalization',
 tooltip=cutoff_tooltip)
max_score = input.float(
 defval=3, title='Max Score Cutoff', group='Physics Normalization',
 tooltip=cutoff_tooltip)

is_auto_norm = input.bool(
 defval=true, title='Auto-Normalization', group='Auto-Normalization',
 tooltip=autonorm_tooltip)
min_auto = input.float(defval=0, title='Min Stdev Multiplier', group='Auto-Normalization')
max_auto = input.float(defval=3, title='Max Stdev Multiplier', group='Auto-Normalization')

gamma = input.float(defval=2.0, title='Gamma', group='Power-Law Normalization', 
 tooltip=gamma_tooltip)

smooth_method = input.string(defval='WMA', title='Method', options=['SMA', 'EMA', 'SMMA (RMA)', 'WMA'],
 group='Smoothing')
smooth_length = input.int(defval=1, title='Length', minval=1, group='Smoothing')

is_plot = input.bool(
 defval=false, title='Plot Window  ', inline='plot0', 
 group='Data Validation')
plot_index = input.int(
 defval=10, title='  ', minval=0, inline='plot0', group='Data Validation')
plot_color = input.color(
 defval=#9598a1, title='', inline='plot0', group='Data Validation', 
 tooltip=plot_dist_tooltip)

// plot_width = input.int(
//  defval=400, title='Width', inline='plot2', group='Data Validation')
// plot_height = input.int(
//  defval=5, title='Height', inline='plot2', group='Data Validation')
// plot_y_shift = input.int(
//  defval=0, title='y-Shift', inline='plot3', group='Data Validation')
// plot_x_shift = input.int(
//  defval=50, title='x-Shift', inline='plot3', group='Data Validation')

plot_width=400
plot_height=5
plot_y_shift=0
plot_x_shift=35

is_print_stats = input.bool(
 defval=false, title='Print Statistics  ', inline='print', group='Data Validation')
is_print_windows = input.bool(
 defval=false, title='Print Windows', inline='print', group='Data Validation')
print_color = input.color(
 defval=color.new(#434651, 0), title='', inline='print', group='Data Validation')

res = input.int(defval=8, title='Resolution', minval=1, group='Multiplicity',
 tooltip=resolution_tooltip)
ind_num = input.int(
 defval=0, title='Stack Number', minval=0, group='Multiplicity',
 tooltip=stack_tooltip)
num_overlays = input.float(
 defval=1, title='Number of Overlays', minval=0, group='Multiplicity',
 tooltip=overlay_tooltip)
shift_percent = input.float(
 defval=0, title='Overlay Shift (%)', group='Multiplicity')

hide_0 = input.bool(
 defval=false, title='Hide', inline='h0', group='Multiplicity')
session_0 = input.session(
 defval="0000-0000", title="", inline='h0', group='Multiplicity',
 tooltip=hide_tooltip)
hide_1 = input.bool(
 defval=false, title='Hide', inline='h1', group='Multiplicity')
session_1 = input.session(
 defval="0000-0000", title="", inline='h1', group='Multiplicity', 
 tooltip=hide_tooltip)
hide_2 = input.bool(
 defval=false, title='Hide', inline='h2', group='Multiplicity')
session_2 = input.session(
 defval="0000-0000", title="", inline='h2', group='Multiplicity',
 tooltip=hide_tooltip)

divide_transp(transp, divisions) =>
    100 * math.pow(transp / 100, 1 / divisions)

filter_color(_color, num_overlays, use_global_coloring, global_color) =>
    color c = _color
    if use_global_coloring
        c := global_color
    transp = color.t(_color)
    if num_overlays > 1
        transp := divide_transp(transp, num_overlays)
    color.new(c, transp)

min_global_color_up = input.color(defval=MIN_COLOR_UP, 
 inline='min', title='Minimum Global Color      Up', group='Global Coloring')
min_global_color_dn = input.color(defval=MIN_COLOR_DN, 
 inline='min', title='Down', group='Global Coloring')
max_global_color_up = input.color(defval=MAX_COLOR_UP, 
 inline='max', title='Maximum Global Color      Up', group='Global Coloring')
max_global_color_dn = input.color(defval=MAX_COLOR_DN, 
 inline='max', title='Down', group='Global Coloring')
use_global_coloring = input.bool(defval=false, 
 title='Use Global Coloring', group='Global Coloring')

min_body_color_up = filter_color(input.color(defval=MIN_COLOR_UP, 
 inline='min0', title='Minimum Body Color       Up', group='Segment Coloring'), 
 num_overlays, use_global_coloring, min_global_color_up)
min_body_color_dn = filter_color(input.color(defval=MIN_COLOR_DN, 
 inline='min0', title='Down', group='Segment Coloring'), num_overlays, 
 use_global_coloring, min_global_color_dn)
max_body_color_up = filter_color(input.color(defval=MAX_COLOR_UP, 
 inline='max0', title='Maximum Body Color       Up', group='Segment Coloring'), 
 num_overlays, use_global_coloring, max_global_color_up)
max_body_color_dn = filter_color(input.color(defval=MAX_COLOR_DN, 
 inline='max0', title='Down', group='Segment Coloring'), num_overlays, 
 use_global_coloring, max_global_color_dn)

min_tail_color_up = filter_color(input.color(defval=MIN_COLOR_UP, 
 inline='min1', title='Minimum Tail Color        Up', group='Segment Coloring'), 
 num_overlays, use_global_coloring, min_global_color_up)
min_tail_color_dn = filter_color(input.color(defval=MIN_COLOR_DN, 
 inline='min1', title='Down', group='Segment Coloring'), num_overlays, 
 use_global_coloring, min_global_color_dn)
max_tail_color_up = filter_color(input.color(defval=MAX_COLOR_UP, 
 inline='max1', title='Maximum Tail Color        Up', group='Segment Coloring'), 
 num_overlays, use_global_coloring, max_global_color_up)
max_tail_color_dn = filter_color(input.color(defval=MAX_COLOR_DN, 
 inline='max1', title='Down', group='Segment Coloring'), 
 num_overlays, use_global_coloring, max_global_color_dn)

min_border_color_up = filter_color(input.color(defval=MIN_COLOR_UP, 
 inline='min2', title='Minimum Border Color     Up', group='Segment Coloring'), 
 num_overlays, use_global_coloring, min_global_color_up)
min_border_color_dn = filter_color(input.color(defval=MIN_COLOR_DN, 
 inline='min2', title='Down', group='Segment Coloring'), num_overlays, 
 use_global_coloring, min_global_color_dn)
max_border_color_up = filter_color(input.color(defval=MAX_COLOR_UP, 
 inline='max2', title='Maximum Border Color     Up', group='Segment Coloring'), 
 num_overlays, use_global_coloring, max_global_color_up)
max_border_color_dn = filter_color(input.color(defval=MAX_COLOR_DN, 
 inline='max2', title='Down', group='Segment Coloring'), num_overlays, 
 use_global_coloring, max_global_color_dn)

min_wick_color_up = filter_color(input.color(defval=MIN_COLOR_UP, 
 inline='min4', title='Minimum Wick Color       Up', group='Segment Coloring'), 
 num_overlays, use_global_coloring, min_global_color_up)
min_wick_color_dn = filter_color(input.color(defval=MIN_COLOR_DN, 
 inline='min4', title='Down', group='Segment Coloring'), num_overlays, 
 use_global_coloring, min_global_color_dn)
max_wick_color_up = filter_color(input.color(defval=MAX_COLOR_UP, 
 inline='max4', title='Maximum Wick Color       Up', group='Segment Coloring'), 
 num_overlays, use_global_coloring, max_global_color_up)
max_wick_color_dn = filter_color(input.color(defval=MAX_COLOR_DN, 
 inline='max4', title='Down', group='Segment Coloring'), num_overlays, 
 use_global_coloring, max_global_color_dn)

min_outline_body_color_up = filter_color(input.color(defval=MIN_COLOR_UP, 
 inline='min3', title='Minimum Border Color     Up', group='Candle Outline Coloring'), 
 num_overlays, use_global_coloring, min_global_color_up)
min_outline_body_color_dn = filter_color(input.color(defval=MIN_COLOR_DN, 
 inline='min3', title='Down', group='Candle Outline Coloring'), num_overlays, 
 use_global_coloring, min_global_color_dn)
max_outline_body_color_up = filter_color(input.color(defval=MAX_COLOR_UP, 
 inline='max3', title='Maximum Border Color     Up', group='Candle Outline Coloring'), 
 num_overlays, use_global_coloring, max_global_color_up)
max_outline_body_color_dn = filter_color(input.color(defval=MAX_COLOR_DN, 
 inline='max3', title='Down', group='Candle Outline Coloring'), num_overlays, 
 use_global_coloring, max_global_color_dn)

min_outline_wick_color_up = filter_color(input.color(defval=MIN_COLOR_UP, 
 inline='min4', title='Minimum Wick Color       Up', group='Candle Outline Coloring'), 
 num_overlays, use_global_coloring, min_global_color_up)
min_outline_wick_color_dn = filter_color(input.color(defval=MIN_COLOR_DN, 
 inline='min4', title='Down', group='Candle Outline Coloring'), num_overlays, 
 use_global_coloring, min_global_color_dn)
max_outline_wick_color_up = filter_color(input.color(defval=MAX_COLOR_UP, 
 inline='max4', title='Maximum Wick Color       Up', group='Candle Outline Coloring'), 
 num_overlays, use_global_coloring, max_global_color_up)
max_outline_wick_color_dn = filter_color(input.color(defval=MAX_COLOR_DN, 
 inline='max4', title='Down', group='Candle Outline Coloring'), num_overlays, 
 use_global_coloring, max_global_color_dn)

close_based = input.bool(defval=false, 
 title='Color by previous close', group='Candle Settings')
fat_wicks = input.bool(defval=false, 
 title='Thick Wicks', group='Candle Settings')

source_to_ticker(source) =>
    pos = str.pos(source, ':')
    prefix = str.substring(source, 0, pos)
    suffix = str.substring(source, pos+1)
    ticker = ticker.new(prefix, suffix, syminfo.session)

chart_ticker= ticker.new(syminfo.prefix, syminfo.ticker, syminfo.session)
price_ticker = source_to_ticker(price_source)
value_ticker = source_to_ticker(physics_source)

[open_, high_, low_, close_] = request.security(
 price_source == '' ? chart_ticker: str.contains(price_source, '-')
 or str.contains(price_source, '+') or str.contains(price_source, '*')
 or str.contains(price_source, '/') ? price_source : price_ticker, 
 timeframe.period, [open, high, low, close])

[iopen, ihigh, ilow, iclose] = request.security_lower_tf(
 price_source == '' ? chart_ticker: str.contains(price_source, '-')
 or str.contains(price_source, '+') or str.contains(price_source, '*')
 or str.contains(price_source, '/') ? price_source : price_ticker, 
 sample_tf, [open, high, low, close])

if timeframe.in_seconds(timeframe.period) >= timeframe.in_seconds("D")
 and array.size(iopen) > 0
    open_ := array.get(iopen, 0)
    high_ := array.max(ihigh)
    low_ := array.min(ilow)
    close_ := array.get(iclose, array.size(iclose)-1)

print(price, _text, _color=color.gray, _textcolor=color.silver, 
 _size=size.auto, _style=label.style_label_down, _yloc=yloc.price) =>
    label.new(bar_index, price, color=_color, 
     textcolor=_textcolor, size=_size, style=_style, 
     yloc=_yloc, text=str.tostring(_text))

get_price_range(bin_num, res) =>
    float max_price = ((bin_num + 1) / res) * (high_ - low_) + low_
    float min_price = (bin_num / res) * (high_ - low_) + low_
    [min_price, max_price]
    
get_body_range() =>
    float body_low = math.min(open_, close_)
    float body_high = math.max(open_, close_)
    [body_low, body_high]

updn_color(zscore, min_value, max_value, min_color_dn, max_color_dn, 
 min_color_up, max_color_up, close_based) =>
    color _color = color.from_gradient(
     zscore, min_value, max_value, min_color_dn, max_color_dn)
    if (close_[0] >= close_[1] and close_based) or (close_ >= open_ and not close_based)
        _color := color.from_gradient(
         zscore, min_value, max_value, min_color_up, max_color_up)
    _color

get_candle_colors(
 zscores, bin_num, res, fat_wicks, close_based, min_value, max_value,
 min_body_color_up, min_body_color_dn, max_body_color_up, max_body_color_dn,
 min_tail_color_up, min_tail_color_dn, max_tail_color_up, max_tail_color_dn) =>
    color _color = EMPTY_COLOR
    [min_price, max_price] = get_price_range(bin_num, res)
    [body_low, body_high] = get_body_range()
    
    float z = 0
    if bin_num < res
        z := array.get(zscores, bin_num)
    
    // if it doesn't have a thin wick
    if not (max_price <= body_low or min_price >= body_high)
        _color := updn_color(z, min_value, max_value, min_body_color_dn, 
         max_body_color_dn, min_body_color_up, max_body_color_up, close_based)
    // if it does have thin wick, but we want a fat tail
    if fat_wicks and (max_price <= body_low or min_price >= body_high)
        _color := updn_color(z, min_value, max_value, min_tail_color_dn, 
         max_tail_color_dn, min_tail_color_up, max_tail_color_up, close_based)
    _color
    
get_wick_color(
 zscores, bin_num, res, fat_wicks, close_based, min_value, max_value,
 min_color_up, min_color_dn, max_color_up, max_color_dn) =>
    color _color = EMPTY_COLOR
    [min_price, max_price] = get_price_range(bin_num, res)
    [body_low, body_high] = get_body_range()
    
    // if it has a thin wick
    if (max_price < body_low or min_price > body_high or 
     (max_price >= body_low and min_price < body_low) or 
     (min_price <= body_high and max_price > body_high)) and not fat_wicks
        float z = 0
        if bin_num < res
            z := array.get(zscores, bin_num)
        _color := updn_color(z, min_value, max_value, min_color_dn, 
         max_color_dn, min_color_up, max_color_up, close_based)
    _color

// match border segment color with body color
get_border_color(
 zscores, bin_num, res, fat_wicks, close_based, min_value, max_value, 
 min_color_up, min_color_dn, max_color_up, max_color_dn) =>
    color _color = EMPTY_COLOR
    [body_low, body_high] = get_body_range()
    [min_price, max_price] = get_price_range(bin_num, res)
    // if inside body
    if not (max_price <= body_low or min_price >= body_high)
        float z = 0
        if bin_num < res
            z := array.get(zscores, bin_num)
        _color := updn_color(z, min_value, max_value, min_color_dn, 
         max_color_dn, min_color_up, max_color_up, close_based)
        _color := color.new(_color, divide_transp(color.t(_color), 2))
    // if outside body but we want a fat wick
    if fat_wicks and (max_price <= body_low or min_price >= body_high)
        _color := EMPTY_COLOR
    _color
    
get_bin_ohlc(bin_num, res, fat_wicks) =>
    [min_price, max_price] = get_price_range(bin_num, res)
    float _open = open_
    float _high = max_price
    float _low = min_price
    float _close = close_

    body_low = math.min(open_, close_)
    body_high = math.max(open_, close_)
    
    // create lower wick
    if max_price < body_low
        _close := max_price
        _open := fat_wicks ? min_price : max_price
    // create lower wick and body
    else if max_price >= body_low and min_price < body_low
     and max_price <= body_high
        _close := max_price
        _open := fat_wicks ? min_price : body_low
    // create body
    else if max_price <= body_high and min_price >= body_low
        _open := min_price
        _close := max_price
    // create upper wick and body
    else if max_price >= body_high and min_price < body_high
     and min_price >= body_low
        _open := min_price
        _close := fat_wicks ? max_price : body_high
    // create upper wick
    else if min_price > body_high
        _open := min_price
        _close := fat_wicks ? max_price : min_price

    [_open, _high, _low, _close]

//http://www.burtonsys.com/climate/composite_standard_deviations.html
update_group_stats(new_count, old_count, new_mean, 
 old_mean, new_stdev, old_stdev) =>
    count = new_count + old_count
    mean = ((new_count * new_mean) + (old_count * old_mean)) / count
    var1 = math.pow(new_stdev, 2)
    essg1 = (new_count - 1) * var1
    var2 = math.pow(old_stdev, 2)
    essg2 = (old_count - 1) * var2
    ess = essg1 + essg2
    gss1 = math.pow(new_mean - mean, 2) * new_count
    gss2 = math.pow(old_mean - mean, 2) * old_count
    tgss = gss1 + gss2
    gv = (ess + tgss) / (count - 1)
    stdev = math.sqrt(gv)
    [count, mean, stdev]

set_value(value, values, index) =>
    prev_value = array.get(values, index)
    if na(prev_value)
        array.set(values, index, value) 
    else
        array.set(values, index, value + prev_value) 

interpolate_segment_values(lows, highs, values, resolution) =>
    float[] segment_values = array.new_float(resolution, na)
    for v=0 to array.size(values)-1
    
        _low = array.get(lows, v)
        _high = array.get(highs, v)
        _val = array.get(values, v)

        lowest = array.min(lows)
        highest = array.max(highs)
        
        seg_low = (_low - low_)/(high_ - low_) * resolution
        seg_high = (_high - low_)/(high_ - low_) * resolution
        
        start = math.floor(seg_low)
        end = math.ceil(seg_high)
        sl = seg_low
        seg_tot = seg_high - seg_low
        
        //
        if na(start) or na(end)
            for i=0 to resolution-1
                set_value(_val / resolution, segment_values, i)

        if start >=0 and end <= resolution
            if (start == end or _low == _high)
                i = int((_low - low_)/(high_ - low_) * resolution)
                i := math.min(i, resolution-1)
                set_value(_val, segment_values, i)    
            else
                for i=start to end-1
                    sh = seg_high
                    if seg_high >= (i+1)
                        sh := i+1
                    seg_val = ((sh - sl) / seg_tot) * _val
                    set_value(seg_val, segment_values, i)      
                    sl := i+1
    segment_values

get_tick_values(value, newbar_value, resolution) =>
    varip float[] highs = array.new_float(0)
    varip float[] lows = array.new_float(0)
    varip float[] closes = array.new_float(0)
    varip float[] values = array.new_float(0)
    varip float prev_value = newbar_value
    
    varip reset = false
    if barstate.isnew
        array.clear(highs)
        array.clear(lows)
        array.clear(closes)
        array.clear(values)
        prev_value := newbar_value
        reset := true

    last = open_
    if array.size(closes) > 0
        last := array.get(closes, array.size(closes)-1)

    if close_ > last
        array.push(highs, close_)
        array.push(lows, last)
    else
        array.push(highs, last)
        array.push(lows, close_)
        
    array.push(closes, close_)
    array.push(values, value - prev_value)

    prev_value := value
    segment_values = interpolate_segment_values(lows, highs, values, resolution)
    segment_values

filter_nans(a) =>
    float[] filter = array.new_float(0)
    for i=0 to array.size(a)-1
        if not na(array.get(a, i))
            array.push(filter, array.get(a, i))
    filter

update_stats(
 is_update, samples, index, prev_index, num_indices, index_change, 
 hist_counts, hist_means, hist_stdevs) =>

    var int window_count = na
    var float window_mean = na
    var float window_stdev = na
    
    _samples = filter_nans(samples)
        

    if (index_change or num_indices == 2) 
     and window_count > 0 and is_update

        hist_count = array.get(hist_counts, prev_index)
        hist_mean = array.get(hist_means, prev_index)
        hist_stdev = array.get(hist_stdevs, prev_index)
    
        if na(hist_count)
            hist_count := window_count
            hist_mean := window_mean
            hist_stdev := window_stdev
        else
            [update_count, update_mean, update_stdev] = update_group_stats(
             window_count, hist_count, window_mean, hist_mean, 
             window_stdev, hist_stdev)
            hist_count := update_count
            hist_mean := update_mean
            hist_stdev := update_stdev

        array.set(hist_counts, prev_index, hist_count)
        array.set(hist_means, prev_index, hist_mean)
        array.set(hist_stdevs, prev_index, hist_stdev)
        
        window_count := na
        window_mean := na
        window_stdev := na

    candle_count = array.size(_samples)
    candle_mean = array.avg(_samples)
    candle_stdev = array.stdev(_samples, biased=false)
    
    if na(candle_stdev)
        candle_stdev := 0.0
    
    if candle_count > 0 and is_update
    
        if na(window_count)
            window_count := candle_count
            window_mean := candle_mean
            window_stdev := candle_stdev
        else
            [update_count, update_mean, update_stdev] = update_group_stats(
             candle_count, window_count, candle_mean, window_mean, 
             candle_stdev, window_stdev)
             
            window_count := update_count
            window_mean := update_mean
            window_stdev := update_stdev

    
calc_zscore(x, mean, stdev) =>
    zscore = (x - mean) / stdev

linears_to_logs(linears) =>
    float[] logs = array.new_float(array.size(linears), 0)
    for i=0 to array.size(linears)-1
        val = array.get(linears, i)
        array.set(logs, i, math.log(val))
    logs

linears_to_powers(linears, gamma) =>
    float[] powers = array.new_float(array.size(linears), 0)
    for i=0 to array.size(linears)-1
        val = array.get(linears, i)
        array.set(powers, i, math.pow(val, gamma))
    powers

linears_to_zscores(linears, mean, stdev) =>
    float[] zscores = array.new_float(array.size(linears), 0)
    for i=0 to array.size(linears)-1
        val = array.get(linears, i)
        float z = calc_zscore(val, mean, stdev)
        array.set(zscores, i, z)
    
    zscores

map_linears(map_type, linears, mean, stdev, gamma) =>
    float[] values = linears
    if map_type == 'Z-Score'
        values := linears_to_zscores(linears, mean, stdev)
    else if map_type == 'Logarithm'
        values := linears_to_logs(linears)
    else if map_type == 'Power-Law'
        values := linears_to_powers(linears, gamma)
    values

extract_timeframe(tf) =>
    unit = str.substring(tf, str.length(tf)-1)
    multiplier = 0
    period = ''
    if unit == 'D' or unit == 'W' or unit == 'M'
        multiplier := 
         int(str.tonumber(str.substring(tf, 0, str.length(tf)-1)))
        period := unit
    [multiplier, period]

verticalize_scores(scores) =>
    string verticals = 'Scores : \n'
    for i=array.size(scores)-1 to 0
        verticals += 
         str.tostring(
          str.format("{0, number, .####}",array.get(scores, i)) + '\n')
    verticals

get_index(_time, cycle_tf, window_tf) =>

    index = -1
    cycle_tf_mins = timeframe.in_seconds(cycle_tf) / 60
    window_tf_mins = timeframe.in_seconds(window_tf) / 60
    
    if cycle_tf_mins <= timeframe.in_seconds("1D") / 60
        index := math.floor(
         (hour(_time) * 60 + minute(_time)) / window_tf_mins)
         
    else if cycle_tf_mins <= (timeframe.in_seconds("1W") / 60)
        index:= math.floor(
         ((dayofweek(_time)-1) * 24 * 60 + hour(_time) * 60 + minute(_time)) 
          / window_tf_mins)

    else if cycle_tf_mins <= timeframe.in_seconds("1M") / 60

        index := math.floor(
         ((dayofmonth(_time)-1) * 24 * 60 + hour(_time) * 60 + minute(_time)) 
          / window_tf_mins)
        
    else if cycle_tf_mins == timeframe.in_seconds("3M") / 60 
     and window_tf_mins % (timeframe.in_seconds("W") / 60) == 0 
     and window_tf_mins < timeframe.in_seconds("M") / 60
        quarter = math.floor((month(_time)-1) / 3) + 1
        first_month = quarter * 3 - 2
        t0_s = timestamp(year(_time), first_month, 0) / 1000
        t1_s = timestamp(year(_time), month(_time), dayofmonth(_time)) / 1000
        index := math.floor((t1_s - t0_s - 1) / timeframe.in_seconds("1W"))
        
    else if cycle_tf_mins <= timeframe.in_seconds("12M") / 60 
     and timeframe.in_seconds(window_tf) % timeframe.in_seconds("M") == 0
        [cyc_mult, _] = extract_timeframe(cycle_tf)
        [win_mult, _] = extract_timeframe(window_tf)
        index := math.floor(((month(_time)-1) % (cyc_mult)) / win_mult)
        
    if index == -1
        runtime.error(
         "The chart and window time frame combination is not available")
    index 

update_mean(count, mean, value) =>
    _count = count + 1
    delta = value - mean
    _mean = delta / _count + mean
    [_count, _mean]

partition_candles(lows, highs, values) =>
    low_highs = array.concat(array.copy(lows), highs)
    array.sort(low_highs)
    parts = array.from(array.get(low_highs, 0))
    // remove duplicates
    for i=1 to array.size(low_highs)-1
        v1 = array.get(low_highs, i)
        v0 = array.get(parts, array.size(parts)-1)
        if v1 != v0
            array.push(parts, v1)
    parts

average_partitions(parts, lows, highs, values) =>

    counts = array.new_float(array.size(parts)-1, 0)
    part_avgs = array.new_float(array.size(parts)-1, 0)

    indices = array.new_int(0)
    for i=0 to array.size(values)-1
        array.push(indices, i)

    part_lows = array.new_float(0)
    part_highs = array.new_float(0)
    for i=0 to array.size(parts)-2
        part_low = array.get(parts, i)
        array.push(part_lows, part_low)
        part_high = array.get(parts, i+1)
        array.push(part_highs, part_high)
        new_indices = array.new_int(0)
        for j=0 to array.size(indices)-1
            k = array.get(indices, j)
            value_low = array.get(lows, k)
            value_high = array.get(highs, k)
            value = array.get(values, k)
            if value_low <= part_low and value_high >= part_high
                [c, m] = update_mean(array.get(counts, i), array.get(part_avgs, i), value)
                array.set(counts, i, c)
                array.set(part_avgs, i, m)
            if value_high > part_high
                array.push(new_indices, k)
        indices := new_indices

    [part_lows, part_highs, part_avgs]

scale_partition_values(part_lows, part_highs, part_values, res) =>
    float segment_size = (high_ - low_) / res
    float[] scaled_parts = array.new_float(0)
    for i=0 to array.size(part_values)-1
        v = array.get(part_values, i)
        h = array.get(part_highs, i)
        l = array.get(part_lows, i)
        array.push(scaled_parts, v * (h-l)/segment_size)
    scaled_parts

percent_change(v0, v1) =>
    percent = 100 * (v1 - v0) / v0

get_velocity(v0, v1) =>
    percent_change(v0, v1)

get_acceleration(v0, v1) =>
    percent_change(v0, v1) - percent_change(v0[1], v1[1])

get_physics(mass_option, physics_option) =>

    float mass = na
    if mass_option == "Volume"
        mass := volume
    else if mass_option == "Time"
        mass := time_close - time

    float motion = na
    motion := switch physics_option
        'Position' => close
        'Velocity' => get_velocity(open, close)
        'Momentum' => get_velocity(open, close)
        'Kinetic Energy' => get_velocity(open, close)
        'Acceleration' => get_acceleration(open, close)
        'Force' => get_acceleration(open, close)
        'Potential Energy' => 1 / get_velocity(open, close)

    [mass, motion, time]

should_have_motion(option) =>
    bool is_sum = switch option
        'Position' => true
        'Velocity' => true
        'Momentum' => true
        'Kinetic Energy' => true
        'Acceleration' => true
        'Force' => true
        'Potential Energy' => true
        => false

should_have_mass(option) =>
    bool is_avg = switch option
        'Mass' => true
        // if mass is said to "have" energy, then time also 
        // has mass via E=mc^2 where time ≈ distance * sqrt(mass/E). 
        'Time' => true
        'Volume' =>  true
        'Momentum' => true
        'Kinetic Energy' => true
        'Force' => true
        'Potential Energy' => true
        => false

multiply_arrays(a0, a1) =>
    float[] product = array.new_float(0)
    for i=0 to array.size(a0)-1
        array.push(product, array.get(a0, i) * array.get(a1, i))
    product

calc_segment_physics(mass_segments, motion_segments, physics_selection, resolution) =>
    float[] segment_values = switch physics_selection
        'Mass' => mass_segments
        'Time' => mass_segments
        'Volume' => mass_segments
        'Position' => motion_segments
        'Velocity' => motion_segments
        'Momentum' => multiply_arrays(mass_segments, motion_segments)
        'Kinetic Energy' => multiply_arrays(mass_segments, 
         multiply_arrays(motion_segments, motion_segments))
        'Acceleration' => motion_segments
        'Force' => multiply_arrays(mass_segments, motion_segments)
        'Potential Energy' => multiply_arrays(mass_segments, motion_segments)
        => array.new_float(resolution, na)

absolute_array(a) =>
    for i=0 to array.size(a)-1
        array.set(a, i, math.abs(array.get(a, i)))
    a

get_segment_values(is_get, mass, masses, motions, lows, highs, 
 sample_tf, resolution, mass_selection, physics_selection) =>

    float[] mass_segments = array.new_float(resolution, na)
    float[] motion_segments = array.new_float(resolution, na)

    if is_get

        has_mass = should_have_mass(physics_selection)
        has_motion = should_have_motion(physics_selection)

        if has_motion and array.size(motions) == array.size(lows)
         and array.size(motions) == array.size(highs)
            partitions = partition_candles(lows, highs, motions)
            if array.size(partitions) > 1
                [part_lows, part_highs, part_avgs] = average_partitions(
                 partitions, lows, highs, motions)
                scaled_values = scale_partition_values(
                 part_lows, part_highs, part_avgs, resolution)
                motion_segments := interpolate_segment_values(
                 part_lows, part_highs, scaled_values, resolution)
            else
                motion_segments := interpolate_segment_values(lows, highs, motions, resolution)
                    
        if has_mass and array.size(masses) == array.size(lows)
         and array.size(masses) == array.size(highs)
            values = masses
            if barstate.isrealtime
                if mass_selection == "Time"
                    mass_segments := get_tick_values(timenow, time, resolution)
                else if mass_selection == "Volume"
                    mass_segments := get_tick_values(mass, 0, resolution)
            else
                mass_segments := interpolate_segment_values(lows, highs, masses, resolution)

    motion_segments := absolute_array(motion_segments)
    segment_values = calc_segment_physics(
     mass_segments, motion_segments, physics_selection, resolution)

get_scores(
 is_get, segment_values, lows, highs, 
 resolution, mean, stdev, norm_type, gamma) =>

    float[] scores = array.new_float(resolution, na)

    if is_get
        scores := map_linears(norm_type, segment_values, mean, stdev, gamma)

    scores

auto_norm(mean, stdev, min_mult, max_mult) =>
    min_cutoff = mean + min_mult*stdev
    max_cutoff = mean + max_mult*stdev
    [min_cutoff, max_cutoff]

print_statistics(
 is_print, scores, index, num_indices, cycle_tf, mean, stdev,
 raw_sample_count, score_sample_count, print_color) =>
    if (ta.change(index) or (timeframe.change(cycle_tf) and num_indices == 2)) 
     and is_print
     
        print(low_, verticalize_scores(scores), _color=print_color, 
         _style=label.style_label_up)
    
        stats = 'Window [' + str.tostring(index) + ']\n'
         + 'Mean:' + '\n' + str.tostring(
         str.format("{0, number, .####}", mean)) 
         + '\n' + 'Std Dev:' + '\n' + str.tostring(
          str.format("{0, number, .####}", stdev)) 
         + '\n' + 'Raw Samples:' + '\n' + str.tostring(
          raw_sample_count)
         + '\n' + 'Score Samples:' + '\n' + str.tostring(
          score_sample_count)
         
        print(high_, stats, _color=print_color)

print_windows(is_print, index, num_indices, cycle_tf, print_color) =>
    if (ta.change(index) or (timeframe.change(cycle_tf) and num_indices == 2)) 
     and is_print
        print(high_, index, _color=print_color, _yloc=yloc.price)

label_sampling_start_date(is_label, print_color) =>
    var bool show = false
    if not show and is_label
        show := true
        print(low_, "Start Data Collection", _color=print_color, _size=size.small)

plot_stat_line(stat_x0, x0, y0, y1, txt, plot_color, label_style) =>
    line li = line.new(stat_x0 + x0, y0, stat_x0 + x0, y1, color=plot_color, 
     style=line.style_dotted, xloc=xloc.bar_index)
    y = y1
    if label_style == label.style_label_down
        y := y0
    label la = label.new(stat_x0 + x0, y, text=txt, size=size.small, 
     color=EMPTY_COLOR, textcolor=plot_color, 
     style=label_style, xloc=xloc.bar_index)
    [li, la]

bin_distribution(samples, min_value, max_value, num_bins) =>
    float[] bins = array.new_float(num_bins, 0)
    bin_size = (max_value - min_value) / num_bins
    for i=0 to array.size(samples)-1
        sample = array.get(samples, i)
        if sample >= min_value and sample <= max_value
            bin_index = math.floor(
             ((sample - min_value) / (max_value - min_value)) * (num_bins-1))
            prev_count = array.get(bins, bin_index)
            array.set(bins, bin_index, prev_count + 1)
    bins

construct_plot(
 x0, y0, x1, y1, stat_x0, stat_x1, samples, plot_index, plot_color, plot_width, 
 text_width, mean, stdev, is_auto_norm, min_auto, max_auto, norm_type, dist_type, 
 min_score, max_score) =>

    var box histogram_frame = na
    var box[] histograms = array.new_box(0)
    var box stat_box = na
    var line[] plot_lines = array.new_line(0)
    var label[] plot_labels = array.new_label(0)

    sample_count = array.size(samples)
    sample_min = array.min(samples)
    sample_max = array.max(samples)
    sample_mean = array.avg(samples)
    sample_stdev = array.stdev(samples, biased=false)
    
    sample_min := sample_min < sample_mean - 4*sample_stdev ? 
     sample_mean - 4*sample_stdev : sample_min
    sample_max := sample_max > sample_mean + 4*sample_stdev ? 
     sample_mean + 4*sample_stdev : sample_max
    
    score_min = min_score
    score_max = max_score
    if is_auto_norm
        [mi, mx] = auto_norm(mean, stdev, min_auto, max_auto)
        score_min := mi
        score_max := mx
    
    bin_width = plot_width - text_width
    
    binned_distribution = bin_distribution(
     samples, sample_min, sample_max, bin_width)

    binned_min = bin_distribution(
     array.from(score_min), sample_min, sample_max, bin_width)
    binned_max = bin_distribution(
     array.from(score_max), sample_min, sample_max, bin_width)
    binned_mean = bin_distribution(
     array.from(sample_mean), sample_min, sample_max, bin_width)
    binned_stdev_pos = bin_distribution(
     array.from(sample_mean+sample_stdev), sample_min, sample_max, bin_width)
    binned_stdev_neg = bin_distribution(
     array.from(sample_mean-sample_stdev), sample_min, sample_max, bin_width)
    
    max_count = array.max(binned_distribution)
    max_height = (y0-y1)

    for i=0 to array.size(binned_distribution)-1
        count = array.get(binned_distribution, i)
        height = (count / max_count) * max_height
    
        array.push(histograms, box.new(
         x0+(i), height + y1, x0+(i+1), y1, xloc=xloc.bar_index, 
         border_color=plot_color, bgcolor=plot_color))
    
    histogram_frame := box.new(x0-1, y0, x1+1, y1, 
     bgcolor=color.new(color.black, 100), 
     text='\n' + norm_type + ' on ' + dist_type + ' Distribution', 
     xloc=xloc.bar_index, text_size=size.normal, text_valign=text.align_top,
     text_halign=text.align_right, border_color=plot_color, text_color=plot_color)
     
    bin_min  = array.indexof(binned_min, 1)
    bin_max = array.lastindexof(binned_max, 1)
    bin_mean = array.indexof(binned_mean, 1)
    bin_stdev_pos = array.indexof(binned_stdev_pos, 1)
    bin_stdev_neg = array.indexof(binned_stdev_neg, 1)

    int[] bin_stats = array.from(
     bin_min, bin_max, bin_mean, bin_stdev_pos, bin_stdev_neg)
    string[] bin_txts = array.from(
     'Min\nCutoff', 'Max\nCutoff', 'Mean', '[+]\nStdev', '[-]\nStdev')
    for i=0 to array.size(bin_stats)-1
        bs = array.get(bin_stats, i)
        bt = array.get(bin_txts, i)
        style = label.style_label_down
        if i > 1
            style := label.style_label_up
            
        line plot_line = na
        label plot_label = na
        if bs != -1
            [li, la] = plot_stat_line(bs, x0, y0, y1, bt, plot_color, style)
            plot_line := li
            plot_label := la
            
        array.push(plot_lines, plot_line)
        array.push(plot_labels, plot_label)

    stat_txt =
     'Window [' + str.tostring(plot_index) + ']' + '\n' +
     'Mean: ' + str.format("{0, Number, .###}", sample_mean) + '\n' + 
     'Std Dev: ' + str.format("{0, Number, .###}", sample_stdev) + '\n' +
     'Min Value: ' + str.format("{0, Number, .###}", sample_min) + '\n' +
     'Max Value: ' + str.format("{0, Number, .###}", sample_max) + '\n' +
     'Min Cutoff: ' + str.format("{0, Number, .###}", score_min) + '\n' +
     'Max Cutoff: ' + str.format("{0, Number, .###}", score_max) + '\n' +
     'Samples: ' + str.tostring(sample_count)
     
    stat_box := box.new(stat_x0, y0, stat_x1, y1, text=stat_txt, 
     bgcolor=EMPTY_COLOR, text_color=plot_color, text_size=size.auto, 
     text_halign=text.align_right, border_color=EMPTY_COLOR)
     
    [histogram_frame, histograms, stat_box, bin_stats, plot_lines, plot_labels]

update_plot(
 x0, x1, stat_x0, stat_x1, histogram_frame, histograms, stat_box, 
 bin_stats, plot_lines, plot_labels) =>

    for i=0 to array.size(histograms)-1
        h = array.get(histograms, i)
        box.set_left(h, x0+i)
        box.set_right(h, x0+(i+1))

    box.set_left(histogram_frame, x0)
    box.set_right(histogram_frame, x1)
    
    box.set_left(stat_box, stat_x0)
    box.set_right(stat_box, stat_x1)
    
    for i=0 to array.size(bin_stats)-1
        bs = array.get(bin_stats, i)
        line.set_xloc(array.get(plot_lines, i), x0+bs, x0+bs, 
         xloc=xloc.bar_index)
        label.set_xloc(array.get(plot_labels, i), x0+bs, 
         xloc=xloc.bar_index)

plot_distribution(
 is_plot, samples, plot_index, plot_color, plot_width, plot_height, 
 plot_x_shift, plot_y_shift, means, stdevs, is_auto_norm, min_auto, 
 max_auto, norm_type, dist_type, min_score, max_score)  =>

    var box histogram_frame = na
    var box[] histograms = na
    var box stat_box = na
    var line[] plot_lines = na
    var label[] plot_labels = na
    var int[] bin_stats = na

    if is_plot
    
        int text_width = 60
        float y0 = close_*(plot_y_shift/100 + 1)
        float y1 = y0 - (plot_height/100)*close_
        int x0 = bar_index + plot_x_shift
        int x1 = x0 + plot_width
        int stat_x0 = x0 + plot_width - text_width
        int stat_x1 = stat_x0+text_width-10

        if barstate.isconfirmed and not na(histograms)
        
            update_plot(
             x0, x1, stat_x0, stat_x1, histogram_frame, histograms, stat_box, 
             bin_stats, plot_lines, plot_labels)
            
        if barstate.islast and na(histograms)
        
            if array.size(samples) == 0
                runtime.error(ERROR_PREFIX + 
                 "No data was collected in this window to plot.  Turn distribution 
                 plotting off or change window number to one that exists in this plane 
                 of existence.")            
            
            mean = array.get(means, plot_index)
            stdev = array.get(stdevs, plot_index)
            
            [hf, hs, sb, bs, lis, las] = 
             construct_plot(
              x0, y0, x1, y1, stat_x0, stat_x1, samples, plot_index, 
              plot_color, plot_width,  text_width, mean, stdev, is_auto_norm, 
              min_auto, max_auto, norm_type, dist_type, min_score, max_score)
        
            histogram_frame := hf
            histograms := hs
            stat_box := sb
            plot_lines := lis
            plot_labels := las
            bin_stats := bs

collect_distribution(is_collect, scores) =>
    var float[] distribution = array.new_float(0)
    if is_collect
        for i=0 to array.size(scores)-1
            if not na(array.get(scores, i))
                array.push(distribution, array.get(scores, i))
    distribution

should_hide_candles(hide_0, session_0, hide_1, session_1, hide_2, session_2) =>
    bool is_hide = false
    if (hide_0 and not na(time(timeframe.period, session_0))) or
     (hide_1 and not na(time(timeframe.period, session_1))) or
     (hide_2 and not na(time(timeframe.period, session_2)))
        is_hide := true
    is_hide

should_collect_samples(use_start_date, start_date, use_end_date, end_date) =>
    var bool is_collect = false
    if (time >= start_date and use_start_date) or not use_start_date
        is_collect := true
    if is_collect and (time >= end_date and use_end_date)
        is_collect := false
    is_collect

zero_array_nans(a) =>
    if array.size(a) > 0
        for i=0 to array.size(a)-1
            array.set(a, i, nz(array.get(a, i)))
    a

calc_sma(v0, v1, length) =>
    delta = v0 - v1
    float sma = delta / length + v1

calc_ema(v0, v1, length) =>
    var alpha = 2 / (length + 1)
    float ema = alpha * v0 + (1 - alpha) * v1

calc_rma(v0, v1, length) =>
    var alpha = 1 / length
    float rma = alpha * v0 + (1 - alpha) * v1

calc_wma(v0, v1, length) =>
    float wma = ((length-1) * v1 + length * v0) / length

smooth_operator(is_smooth, values, value_mean, smooth_method, smooth_length) =>

    smoothed_values = values
    if is_smooth
        rolling_mean = switch smooth_method
            "SMA" => ta.sma(value_mean, smooth_length-1)
            "EMA" => ta.ema(value_mean, smooth_length-1)
            "SMMA (RMA)" => ta.rma(value_mean, smooth_length-1)
            "WMA" => ta.wma(value_mean, smooth_length-1)
        smoothed_values := array.new_float(0)
        for i=0 to array.size(values)-1
            switch smooth_method
                "SMA" => array.push(smoothed_values, 
                 calc_sma(array.get(values, i), rolling_mean, smooth_length))
                "EMA" => array.push(smoothed_values, 
                 calc_ema(array.get(values, i), rolling_mean, smooth_length))
                "SMMA (RMA)" => array.push(smoothed_values, 
                 calc_rma(array.get(values, i), rolling_mean, smooth_length))
                "WMA" => array.push(smoothed_values, 
                 calc_wma(array.get(values, i), rolling_mean, smooth_length))

    smoothed_values

input_validation(sample_tf, window_tf, chart_tf, cycle_tf, use_start_date, 
 start_date, use_end_date, end_date, masses, physics_selection, physics_source) =>
    
    sample_tf_mins = timeframe.in_seconds(sample_tf) / 60
    window_tf_mins = timeframe.in_seconds(window_tf) / 60
    chart_tf_mins = timeframe.in_seconds(timeframe.period) / 60
    cycle_tf_mins = timeframe.in_seconds(cycle_tf) / 60
    
    if window_tf_mins > cycle_tf_mins
        runtime.error(ERROR_PREFIX + 
         "Sample time frame must be lower or equal to the cycle timeframe")
    else if cycle_tf_mins == timeframe.in_seconds("M") / 60 
     and window_tf_mins < timeframe.in_seconds("D") / 60
        runtime.error(ERROR_PREFIX + 
         "Sample time frame is too low for the cycle timeframe")
    else if cycle_tf_mins == timeframe.in_seconds("3M") / 60 
     and window_tf_mins < timeframe.in_seconds("W") / 60
        runtime.error(ERROR_PREFIX + 
         "Sample time frame is too low for the cycle timeframe")
    else if cycle_tf_mins == timeframe.in_seconds("12M") / 60 
     and window_tf_mins < timeframe.in_seconds("M") / 60
        runtime.error(ERROR_PREFIX + 
         "Sample time frame is too low for the cycle timeframe")
    
    if chart_tf_mins < sample_tf_mins
        runtime.error(ERROR_PREFIX + 
         "Chart time frame is too low for sample time frame")
    if chart_tf_mins % sample_tf_mins != 0
        runtime.error(ERROR_PREFIX + 
         "Chart time frame should be evenly divisible by sample time frame")

    if use_start_date and start_date > timenow
        runtime.error(ERROR_PREFIX +
         "The Data Collection Start Date is in the future.  
         Please taxi a time machine to start collecting that data")
    if use_start_date and use_end_date and end_date <= start_date
        runtime.error(ERROR_PREFIX +
         "The Data Collection Start Date is at or after the End Date. 
         Did Bill and Ted take you for a ride again?")

    var float total_mass = 0
    var bool has_mass = should_have_mass(physics_selection)
    if has_mass
        if array.size(masses) > 0
            total_mass += array.sum(zero_array_nans(masses))
        if barstate.islast and (total_mass == 0 or na(total_mass))
            source = physics_source == '' ? syminfo.ticker : physics_source
            runtime.error("No volume is available from " + 
              + source + " to calculate " + physics_selection)

update_minmax(values, mins, maxs, index) =>
    min = array.get(mins, index)
    max = array.get(maxs, index)
    if na(min)
        min := array.get(values, 0)
    if na(max)
        max := array.get(values, 0)
    for i=0 to array.size(values)-1
        v = array.get(values, i)
        if v < min
            min := v
        if v > max
            max := v
    array.set(mins, index, min)
    array.set(maxs, index, max)

sample_tf_mins = timeframe.in_seconds(sample_tf) / 60
window_tf_mins = timeframe.in_seconds(window_tf) / 60
chart_tf_mins = timeframe.in_seconds(timeframe.period) / 60
cycle_tf_mins = timeframe.in_seconds(cycle_tf) / 60
scale_factor = (chart_tf_mins / sample_tf_mins)

num_indices = math.ceil((cycle_tf_mins / window_tf_mins)) + 1

var int[] sample_counts = array.new_int(num_indices, 0)
var float[] seg_counts = array.new_float(num_indices, na)
var float[] seg_means = array.new_float(num_indices, na)
var float[] seg_stdevs = array.new_float(num_indices, na)
var float[] score_counts = array.new_float(num_indices, na)
var float[] score_means = array.new_float(num_indices, na)
var float[] score_stdevs = array.new_float(num_indices, na)
var float[] score_min = array.new_float(num_indices, na)
var float[] score_max = array.new_float(num_indices, na)

[masses, motions, mtimes] = request.security_lower_tf(
 physics_source == '' ? chart_ticker: str.contains(physics_source, '-') 
 or str.contains(physics_source, '+') or str.contains(physics_source, '*') 
 or str.contains(physics_source, '/') ? physics_source : value_ticker , 
 sample_tf, get_physics(mass_selection, physics_selection))

// Ironically, the sum of lower timeframe volumes from security_lower_tf
// is not equal ot the real-time volume so a separate request is needed
[mass, motion, _] = request.security(
 physics_source == '' ? chart_ticker: str.contains(physics_source, '-') 
 or str.contains(physics_source, '+') or str.contains(physics_source, '*') 
 or str.contains(physics_source, '/') ? physics_source : value_ticker , 
 timeframe.period, get_physics(mass_selection, physics_selection))

input_validation(sample_tf, window_tf, timeframe.period, cycle_tf, 
 use_start_date, start_date, use_end_date, end_date, masses, 
 physics_selection, physics_source)

[highs, lows, hltimes] = request.security_lower_tf(
 price_source == '' ? chart_ticker: str.contains(price_source, '-')
 or str.contains(price_source, '+') or str.contains(price_source, '*')
 or str.contains(price_source, '/') ? price_source : price_ticker, 
 sample_tf, [high, low, time])

is_hide = should_hide_candles(
 hide_0, session_0, hide_1, session_1, hide_2, session_2)
 
is_collect = should_collect_samples(
 use_start_date, start_date, use_end_date, end_date)

first_index = get_index(time, cycle_tf, window_tf)
last_index = get_index(time_close-1, cycle_tf, window_tf)

normalize(value, min, max) =>
    v = value < min ? min : value
    v := value > max ? max : v
    float norm = (v - min) / (max - min)
    norm

var prev_index = 0
float[] scores = array.new_float(res, na)
int[] counts = array.new_int(res, 0)
var float plot_min_score = 0
var float plot_max_score = 0
var float[] distribution = na

for index=first_index to last_index

    float[] index_masses = array.new_float(0)
    float[] index_motions = array.new_float(0)
    float[] index_highs = array.new_float(0)
    float[] index_lows = array.new_float(0)

    if array.size(masses) > 0 and array.size(masses) == array.size(lows)

        for i=0 to array.size(masses)-1
            if index == get_index(array.get(mtimes, i), cycle_tf, window_tf) and
             index == get_index(array.get(hltimes, i), cycle_tf, window_tf)
                array.push(index_masses, array.get(masses, i))
                array.push(index_motions, array.get(motions, i))
                array.push(index_highs, array.get(highs, i))
                array.push(index_lows, array.get(lows, i))

        array.set(sample_counts, index, 
         array.get(sample_counts, index) + array.size(index_masses))

        segment_values = get_segment_values(
         array.size(index_masses) > 0, 
         mass, index_masses, index_motions, index_lows, index_highs, sample_tf, 
         res, mass_selection, physics_selection)

        if dist_type == 'Log-Normal'
            segment_values := linears_to_logs(segment_values)

        index_change = false
        if index != prev_index
            index_change := true

        update_stats(
         is_collect, 
         segment_values, index, prev_index, num_indices, index_change,
         seg_counts, seg_means, seg_stdevs)
        
        index_scores = get_scores(
         not is_hide and array.size(index_masses) > 0, 
         segment_values, index_lows, index_highs, res, array.get(seg_means, index), 
         array.get(seg_stdevs, index), norm_type, gamma)

        distribution := collect_distribution(
         is_plot and is_collect and array.size(masses) > 0 and index == plot_index, 
         index_scores)        

        update_stats(
         not is_hide, 
         index_scores, index, prev_index, num_indices, index_change,
         score_counts, score_means, score_stdevs)

        if is_auto_norm
            [mi, mx] = auto_norm(
             array.get(score_means, index), array.get(score_stdevs, index), 
             min_auto, max_auto)
            min_score := mi
            max_score := mx

        if index == plot_index
            plot_min_score := min_score
            plot_max_score := max_score

        for i=0 to array.size(scores)-1
            index_score = array.get(index_scores, i)
            score = array.get(scores, i)
            if not na(index_score) and na(score)
                norm = normalize(index_score, min_score, max_score)
                array.set(counts, i, 1)
                array.set(scores, i, norm)
            else if not na(index_score) and not na(score)
                norm = normalize(index_score, min_score, max_score)
                [c, m] = update_mean(array.get(counts, i), score, norm)
                array.set(counts, i, c)
                array.set(scores, i, m)

         
    prev_index := index

score_mean = array.avg(scores)

for i=0 to array.size(scores)-1
    if na(array.get(scores, i))
        array.set(scores, i, 0)

print_statistics(
 is_print_stats and not is_hide, 
 scores, first_index, num_indices, cycle_tf,
 array.get(score_means, first_index), array.get(score_stdevs, first_index), 
 array.get(sample_counts, first_index), array.get(score_counts, first_index),
 print_color)

print_windows(
 is_print_windows and not is_hide, 
 first_index, num_indices, cycle_tf, print_color)

label_sampling_start_date(
 array.sum(masses) and is_collect, print_color)

plot_distribution(
 is_plot, 
 distribution, plot_index, plot_color, plot_width, plot_height, plot_x_shift, 
 plot_y_shift, score_means, score_stdevs, is_auto_norm, min_auto, max_auto, norm_type, 
 dist_type, plot_min_score, plot_max_score)

scores := smooth_operator(
 smooth_length > 1 and not is_hide and array.size(masses) > 0 and not na(score_mean[1]),
 scores, score_mean[1], smooth_method, smooth_length)

if smooth_length > 1 and not is_hide and array.size(masses) > 0
    score_mean := array.avg(scores)

offset = ind_num * 8
shift = (shift_percent/100)*close

[o0, h0, l0, c0] = get_bin_ohlc(0+offset, res, fat_wicks)
plotcandle(is_hide or res < 1 ? na : o0+shift, h0+shift, l0+shift, c0+shift, 
 color=get_candle_colors(scores, 0+offset, res, fat_wicks, close_based, 
 0, 1, min_body_color_up, min_body_color_dn, max_body_color_up, 
 max_body_color_dn, min_tail_color_up, min_tail_color_dn, max_tail_color_up, 
 max_tail_color_dn),
 wickcolor=get_wick_color(scores, 0+offset, res, fat_wicks, close_based, 
 0, 1, min_wick_color_up, min_wick_color_dn, max_wick_color_up, 
 max_wick_color_dn),
 bordercolor=get_border_color(scores, 0+offset, res, fat_wicks, close_based, 
 0, 1, min_border_color_up, min_border_color_dn,  
 max_border_color_up, max_border_color_dn),
 display=display.all - display.price_scale, editable=false)
[o1, h1, l1, c1] = get_bin_ohlc(1+offset, res, fat_wicks)
plotcandle(is_hide or res < 2 ? na : o1+shift, h1+shift, l1+shift, c1+shift, 
 color=get_candle_colors(scores, 1+offset, res, fat_wicks, close_based, 
 0, 1, min_body_color_up, min_body_color_dn, max_body_color_up, 
 max_body_color_dn, min_tail_color_up, min_tail_color_dn, max_tail_color_up, 
 max_tail_color_dn),
 wickcolor=get_wick_color(scores, 1+offset, res, fat_wicks, close_based, 
 0, 1, min_wick_color_up, min_wick_color_dn, max_wick_color_up, 
 max_wick_color_dn),
 bordercolor=get_border_color(scores, 1+offset, res, fat_wicks, close_based, 
 0, 1, min_border_color_up, min_border_color_dn,  
 max_border_color_up, max_border_color_dn),
 display=display.all - display.price_scale, editable=false)
[o2, h2, l2, c2] = get_bin_ohlc(2+offset, res, fat_wicks)
plotcandle(is_hide or res < 3 ? na : o2+shift, h2+shift, l2+shift, c2+shift, 
 color=get_candle_colors(scores, 2+offset, res, fat_wicks, close_based, 
 0, 1, min_body_color_up, min_body_color_dn, max_body_color_up, 
 max_body_color_dn, min_tail_color_up, min_tail_color_dn, max_tail_color_up, 
 max_tail_color_dn),
 wickcolor=get_wick_color(scores, 2+offset, res, fat_wicks, close_based, 
 0, 1, min_wick_color_up, min_wick_color_dn, max_wick_color_up, 
 max_wick_color_dn),
 bordercolor=get_border_color(scores, 2+offset, res, fat_wicks, close_based, 
 0, 1, min_border_color_up, min_border_color_dn,  
 max_border_color_up, max_border_color_dn),
 display=display.all - display.price_scale, editable=false)
[o3, h3, l3, c3] = get_bin_ohlc(3+offset, res, fat_wicks)
plotcandle(is_hide or res < 4 ? na : o3+shift, h3+shift, l3+shift, c3+shift, 
 color=get_candle_colors(scores, 3+offset, res, fat_wicks, close_based, 
 0, 1, min_body_color_up, min_body_color_dn, max_body_color_up, 
 max_body_color_dn, min_tail_color_up, min_tail_color_dn, max_tail_color_up, 
 max_tail_color_dn),
 wickcolor=get_wick_color(scores, 3+offset, res, fat_wicks, close_based, 
 0, 1, min_wick_color_up, min_wick_color_dn, max_wick_color_up, 
 max_wick_color_dn),
 bordercolor=get_border_color(scores, 3+offset, res, fat_wicks, close_based, 
 0, 1, min_border_color_up, min_border_color_dn,  
 max_border_color_up, max_border_color_dn),
 display=display.all - display.price_scale, editable=false)
[o4, h4, l4, c4] = get_bin_ohlc(4+offset, res, fat_wicks)
plotcandle(is_hide or res < 5 ? na : o4+shift, h4+shift, l4+shift, c4+shift, 
 color=get_candle_colors(scores, 4+offset, res, fat_wicks, close_based, 
 0, 1, min_body_color_up, min_body_color_dn, max_body_color_up, 
 max_body_color_dn, min_tail_color_up, min_tail_color_dn, max_tail_color_up, 
 max_tail_color_dn),
 wickcolor=get_wick_color(scores, 4+offset, res, fat_wicks, close_based, 
 0, 1, min_wick_color_up, min_wick_color_dn, max_wick_color_up, 
 max_wick_color_dn),
 bordercolor=get_border_color(scores, 4+offset, res, fat_wicks, close_based, 
 0, 1, min_border_color_up, min_border_color_dn,  
 max_border_color_up, max_border_color_dn),
 display=display.all - display.price_scale, editable=false)
[o5, h5, l5, c5] = get_bin_ohlc(5+offset, res, fat_wicks)
plotcandle(is_hide or res < 6 ? na : o5+shift, h5+shift, l5+shift, c5+shift, 
 color=get_candle_colors(scores, 5+offset, res, fat_wicks, close_based, 
 0, 1, min_body_color_up, min_body_color_dn, max_body_color_up, 
 max_body_color_dn, min_tail_color_up, min_tail_color_dn, max_tail_color_up, 
 max_tail_color_dn),
 wickcolor=get_wick_color(scores, 5+offset, res, fat_wicks, close_based, 
 0, 1, min_wick_color_up, min_wick_color_dn, max_wick_color_up, 
 max_wick_color_dn),
 bordercolor=get_border_color(scores, 5+offset, res, fat_wicks, close_based, 
 0, 1, min_border_color_up, min_border_color_dn,  
 max_border_color_up, max_border_color_dn),
 display=display.all - display.price_scale, editable=false)
[o6, h6, l6, c6] = get_bin_ohlc(6+offset, res, fat_wicks)
plotcandle(is_hide or res < 7 ? na : o6+shift, h6+shift, l6+shift, c6+shift, 
 color=get_candle_colors(scores, 6+offset, res, fat_wicks, close_based, 
 0, 1, min_body_color_up, min_body_color_dn, max_body_color_up, 
 max_body_color_dn, min_tail_color_up, min_tail_color_dn, max_tail_color_up, 
 max_tail_color_dn),
 wickcolor=get_wick_color(scores, 6+offset, res, fat_wicks, close_based, 
 0, 1, min_wick_color_up, min_wick_color_dn, max_wick_color_up, 
 max_wick_color_dn),
 bordercolor=get_border_color(scores, 6+offset, res, fat_wicks, close_based, 
 0, 1, min_border_color_up, min_border_color_dn,  
 max_border_color_up, max_border_color_dn),
 display=display.all - display.price_scale, editable=false)
[o7, h7, l7, c7] = get_bin_ohlc(7+offset, res, fat_wicks)
plotcandle(is_hide or res < 8 ? na : o7+shift, h7+shift, l7+shift, c7+shift, 
 color=get_candle_colors(scores, 7+offset, res, fat_wicks, close_based, 
 0, 1, min_body_color_up, min_body_color_dn, max_body_color_up, 
 max_body_color_dn, min_tail_color_up, min_tail_color_dn, max_tail_color_up, 
 max_tail_color_dn),
 wickcolor=get_wick_color(scores, 7+offset, res, fat_wicks, close_based, 
 0, 1, min_wick_color_up, min_wick_color_dn, max_wick_color_up, 
 max_wick_color_dn),
 bordercolor=get_border_color(scores, 7+offset, res, fat_wicks, close_based, 
 0, 1, min_border_color_up, min_border_color_dn,  
 max_border_color_up, max_border_color_dn),
 display=display.all - display.price_scale, editable=false)
 
// skeletor candle
plotcandle(is_hide ? na : 
 open_+shift, high_+shift, low_+shift, close_+shift, color=EMPTY_COLOR, 
 wickcolor=ind_num == 0 ? updn_color(score_mean, 0, 1, 
 min_outline_wick_color_dn, max_outline_wick_color_dn, min_outline_wick_color_up, 
 max_outline_wick_color_up, close_based): EMPTY_COLOR, 
 bordercolor=ind_num == 0 ? updn_color(score_mean, 0, 1, 
 min_outline_body_color_dn, max_outline_body_color_dn, min_outline_body_color_up, 
 max_outline_body_color_up, close_based) : EMPTY_COLOR, 
 display=display.all - display.price_scale, editable=false)
 