def plot_robots_time_micmac(deploy_robots_mic, deploy_robots_mac, species_ind, node_ind):
    
    fig = plt.figure()
    #plt.axis('equal')
    num_nodes = deploy_robots_mic.shape[0]
    num_iter = deploy_robots_mic.shape[1]

    delta_t = 0.04
    x = np.arange(0, num_iter) * delta_t


    # plot evolution of robot population over time
    labeled = False
    for n in node_ind:    
        y = deploy_robots_mic[n,:,species_ind]
        y2 = deploy_robots_mac[n,:,species_ind]

        if (not labeled):
            labeled = True
            plt.plot(x,y,color='green', label='Microscopic')
            plt.plot(x,y2, color='red', label='Macroscopic')
        else:    
            plt.plot(x,y,color='green')
            plt.plot(x,y2, color='red')
            
    # plot legend and labels
    plt.legend(loc='upper right', shadow=False, fontsize='large')     
    plt.xlabel('Time [s]')    
    plt.ylabel('Number of robots')
    
    #plt.show()
    return fig

def plotResults(datasetName, sampleSizes, foldsSet, cvScalings, sampleMethods, fileNameSuffix):
    """
    Plots the errors for a particular dataset on a bar graph. 
    """

    for k in range(len(sampleMethods)):
        outfileName = outputDir + datasetName + sampleMethods[k] + fileNameSuffix + ".npz"
        data = numpy.load(outfileName)

        errors = data["arr_0"]
        meanMeasures = numpy.mean(errors, 0)

        for i in range(sampleSizes.shape[0]):
            plt.figure(k*len(sampleMethods) + i)
            plt.title("n="+str(sampleSizes[i]) + " " + sampleMethods[k])

            for j in range(errors.shape[3]):
                plt.plot(foldsSet, meanMeasures[i, :, j])
                plt.xlabel("Folds")
                plt.ylabel('Error')

            labels = ["VFCV", "PenVF+"]
            labels.extend(["VFP s=" + str(x) for x in cvScalings])
            plt.legend(tuple(labels))
    plt.show()

def scree_plot(pca_obj, fname=None): 
    '''
    Scree plot for variance & cumulative variance by component from PCA. 

    Arguments: 
        - pca_obj: a fitted sklearn PCA instance
        - fname: path to write plot to file

    Output: 
        - scree plot 
    '''   
    components = pca_obj.n_components_ 
    variance = pca.explained_variance_ratio_
    plt.figure()
    plt.plot(np.arange(1, components + 1), np.cumsum(variance), label='Cumulative Variance')
    plt.plot(np.arange(1, components + 1), variance, label='Variance')
    plt.xlim([0.8, components]); plt.ylim([0.0, 1.01])
    plt.xlabel('No. Components', labelpad=11); plt.ylabel('Variance Explained', labelpad=11)
    plt.legend(loc='best') 
    plt.tight_layout() 
    if fname is not None:
        plt.savefig(fname)
        plt.close() 
    else:
        plt.show() 
    return 

def plot_scenario(strategies, names, scenario_id=1):
    probabilities = get_scenario(scenario_id)

    plt.figure(figsize=(6, 4.5))

    ax = plt.subplot(111)
    ax.spines["top"].set_visible(False)
    ax.spines["bottom"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.spines["left"].set_visible(False)

    ax.get_xaxis().tick_bottom()
    ax.get_yaxis().tick_left()

    plt.yticks(fontsize=14)
    plt.xticks(fontsize=14)
    plt.xlim((0, 1300))

    # Remove the tick marks; they are unnecessary with the tick lines we just plotted.
    plt.tick_params(axis="both", which="both", bottom="on", top="off",
                    labelbottom="on", left="off", right="off", labelleft="on")

    for rank, (strategy, name) in enumerate(zip(strategies, names)):
        plot_strategy(probabilities, strategy, name, rank)

    plt.title("Bandits: " + str(probabilities), fontweight='bold')
    plt.xlabel('Number of Trials', fontsize=14)
    plt.ylabel('Cumulative Regret', fontsize=14)
    plt.legend(names)
    plt.show()

 def default_run(self):
     """
     Plots the results, saves the figure, and finally displays it from simulating codewords with Sum-prod and Max-prod
     algorithms across variance levels. This combines the results in one plot.
     :return:
     """
     if not os.path.exists("./graphs"):
         os.makedirs("./graphs")
     self.save_time = str(int(time.time()))
     self.simulate(Decoder.SUM_PROD)
     self.compute_error()
     plt.plot([math.log10(x) for x in self.variance_levels], [math.log10(y) for y in self.bit_error_probability],
              "ro-", label="Sum-Prod")
     self.simulate(Decoder.MAX_PROD)
     self.compute_error()
     plt.plot([math.log10(x) for x in self.variance_levels], [math.log10(y) for y in self.bit_error_probability],
              "g^--", label="Max-Prod")
     plt.legend(loc=2)
     plt.title("Hamming Decoder Factor Graph Simulation Results\n" +
               r"$\log_{10}(\sigma^2)$ vs. $\log_{10}(P_e)$" + " for Max-Prod & Sum-Prod Algorithms\n" +
               "Sample Size n = %(codewords)s Codewords \n Variance Levels = %(levels)s"
               % {"codewords": str(self.iterations), "levels": str(self.variance_levels)})
     plt.xlabel("$\log_{10}(\sigma^2)$")
     plt.ylabel(r"$\log_{10}(P_e)$")
     plt.savefig("graphs/%(time)s-max-prod-sum-prod-%(num_codewords)s-codewords-variance-bit_error_probability.png" %
                 {"time": self.save_time,
                  "num_codewords": str(self.iterations)}, bbox_inches="tight")
     plt.show()

  def test_get_obs(self):

    plt.figure()
    ant_sigs = antennas.antennas_signal(self.ants, self.ant_models, self.sources, self.rad.timebase)
    rad_sig_full = self.rad.sampled_signal(ant_sigs[0, :], 0)
    obs_full = self.rad.get_full_obs(ant_sigs, self.utc_date, self.config)

    ant_sigs_simp = antennas.antennas_simplified_signal(self.ants, self.ant_models, self.sources, self.rad.baseband_timebase, self.rad.int_freq)
    obs_simp = self.rad.get_simplified_obs(ant_sigs_simp, self.utc_date, self.config)


    freqs, spec_full_before_obs = spectrum.plotSpectrum(rad_sig_full, self.rad.ref_freq, label='full_before_obs_obj', c='blue')
    freqs, spec_full = spectrum.plotSpectrum(obs_full.get_antenna(1), self.rad.ref_freq, label='full', c='cyan')
    freqs, spec_simp = spectrum.plotSpectrum(obs_simp.get_antenna(1), self.rad.ref_freq, label='simp', c='red')
    plt.legend()

    self.assertTrue((spec_full_before_obs == spec_full).all(), True)


    plt.figure()
    plt.plot(freqs, (spec_simp-spec_full)/spec_full)
    plt.show()

    print len(obs_full.get_antenna(1)), obs_full.get_antenna(1).mean()
    print len(obs_simp.get_antenna(1)), obs_simp.get_antenna(1).mean()

def plotErrorBars(dict_to_plot, x_lim, y_lim, xlabel, y_label, title, out_file, margin=[0.05, 0.05], loc=2):

    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(y_label)

    if y_lim is None:
        y_lim = [1 * float("Inf"), -1 * float("Inf")]

    max_val_seen_y = y_lim[1] - margin[1]
    min_val_seen_y = y_lim[0] + margin[1]
    print min_val_seen_y, max_val_seen_y
    max_val_seen_x = x_lim[1] - margin[0]
    min_val_seen_x = x_lim[0] + margin[0]
    handles = []
    for k in dict_to_plot:
        means, stds, x_vals = dict_to_plot[k]

        min_val_seen_y = min(min(np.array(means) - np.array(stds)), min_val_seen_y)
        max_val_seen_y = max(max(np.array(means) + np.array(stds)), max_val_seen_y)

        min_val_seen_x = min(min(x_vals), min_val_seen_x)
        max_val_seen_x = max(max(x_vals), max_val_seen_x)

        handle = plt.errorbar(x_vals, means, yerr=stds)
        handles.append(handle)
        print max_val_seen_y
    plt.xlim([min_val_seen_x - margin[0], max_val_seen_x + margin[0]])
    plt.ylim([min_val_seen_y - margin[1], max_val_seen_y + margin[1]])
    plt.legend(handles, dict_to_plot.keys(), loc=loc)
    plt.savefig(out_file)

def make_overview_plot(filename, title, noip_arrs, ip_arrs):
    plt.title("Inner parallelism - " + title)

    
    plt.ylabel('Time (ms)', fontsize=12)

    x = 0
    barwidth = 0.5
    bargroupspacing = 1.5

    for z in zip(noip_arrs, ip_arrs):
        noip,ip = z
        noip_mean,noip_conf = conf_stats(noip)
        ip_mean,ip_conf = conf_stats(ip)

        b_noip = plt.bar(x, noip_mean, barwidth, color='r', yerr=noip_conf, ecolor='black', alpha=0.7)
        x += barwidth

        b_ip = plt.bar(x, ip_mean, barwidth, color='b', yerr=ip_conf, ecolor='black', alpha=0.7)
        x += bargroupspacing

    plt.xticks([0.5, 2.5, 4.5], ['50k', '100k', '200k'], rotation='horizontal')

    fontP = FontProperties()
    fontP.set_size('small')

    plt.legend([b_noip, b_ip], \
        ('no inner parallelism', 'inner parallelism'), \
        prop=fontP, loc='upper center', bbox_to_anchor=(0.5, -0.05), fancybox=True, shadow=True, ncol=2)
   
    plt.ylim([0,62000])
    plt.savefig(output_file(filename))
    plt.clf()

def main():
    parser = argparse.ArgumentParser(description="""Compute subset of users who rated at
                                     least 10 movies and plot fraction of users satisfied
                                     as a function of inventory size.""")
    parser.add_argument("infilename",
                        help="Read from this file.", type=open)
    args = parser.parse_args()

    ratings = read_inputs(args.infilename)
    ratings = ratings.drop("timestamp", axis=1)
    movie_rankings = find_movie_rankings(ratings)
    ratings = ratings.drop("rating", axis=1)
    user_rankings = find_user_rankings(ratings, movie_rankings)
    num_users = user_rankings.user_id.unique().size
    num_movies = movie_rankings.shape[0]
    user_rankings = clean_rankings(user_rankings)

    us_levels_100 = find_satisfaction(user_rankings, num_users, num_movies)
    us_levels_90 = find_satisfaction(user_rankings, num_users, num_movies, satisfaction_level=0.9)

    rc('text', usetex=True)
    plt.title('Percent of Users Satisfied vs Inventory Size in the MovieLens Dataset')
    plt.xlabel('Inventory Size')
    plt.ylabel('Percent of Users Satisfied')
    plt.plot(us_levels_100, 'b', label=r'$100\% \ satisfaction$')
    plt.plot(us_levels_90, 'r--', label=r'$90\% \ satisfaction$')
    plt.legend()
    d = datetime.datetime.now().isoformat()
    plt.savefig('user_satisfaction_%s.png' % d)

def create_lib_complexity_plot(fastqFile,mer_size=20):
    
    #def calculate_library_complexity(fastqfiles):
    
    report_uniqness_after_N_reads = 100000
    #middle_mer_start 
    
    randm_mer_fastq_complexity = {}
    starting_mer_fastq_complexity = {}
    count_uniq_random_mers =0
    count_uniq_starting_mers = 0
    
    bucket_percent_uniq_random_mers = []
    bucket_percent_uniq_starting_mers = []
    
    
    #output file name
    #fastq_prefix = get_guessed_fastq_prefix(fastqFile)
    fastq_prefix = fastqFile + '_kmer' + str(mer_size)
    baseDir = os.path.dirname(fastqFile) or '.'
    
    lib_complexity_plot_name = baseDir + '/' + fastq_prefix + '_lib_complexity_plot.png'
    lib_complexity_numbers_file = baseDir + '/' + fastq_prefix + '_lib_complexity.data'

    fh = open(lib_complexity_numbers_file,'w')

    for count,read in enumerate(get_genericfastq_iterator(fastqFile)):
        if count % report_uniqness_after_N_reads == 0 and count > 0:
            uniq_random_mer_percent = round( (count_uniq_random_mers/float(count+1))*100, 2)
            uniq_starting_mer_percent = round( (count_uniq_starting_mers/float(count+1))*100, 2)
            bucket_percent_uniq_random_mers.append(uniq_random_mer_percent)
            bucket_percent_uniq_starting_mers.append( uniq_starting_mer_percent  )
            outline = '%s \t %s \t %s \t %s \t %s \n' %  (count, count_uniq_random_mers,
                                                        count_uniq_starting_mers,
                                                        uniq_random_mer_percent,
                                                        uniq_starting_mer_percent)
            fh.write(outline)
            
        starting_mer = read.seq[:mer_size]
        random_mer_loc=random.randint(0,(len(read.seq)-mer_size))
        random_mer = read.seq[random_mer_loc:random_mer_loc+mer_size]
        
        if randm_mer_fastq_complexity.get(random_mer,None) is None:
            count_uniq_random_mers +=1
            randm_mer_fastq_complexity[random_mer] = 1
        
        if starting_mer_fastq_complexity.get(starting_mer,None) is None:
            count_uniq_starting_mers += 1
            starting_mer_fastq_complexity[starting_mer] = 1
        
    fh.close()
    
    #plotting the library complexity
    x_axis = [ x*report_uniqness_after_N_reads for x in xrange(1,len(bucket_percent_uniq_random_mers)+1) ]
    plt.plot(x_axis,bucket_percent_uniq_random_mers,label='rdm_%dmer' % mer_size)
    plt.plot(x_axis,bucket_percent_uniq_starting_mers,label='start_%dmer' % mer_size)
    plt.ylabel('percent unqiue')
    plt.xlabel('Read Counts')
    plt.legend()
    plt.savefig(lib_complexity_plot_name)

def plot_multiple_likelhood_values(likelihood_arr, time_axis=0, 
                                   x=None, save_path='',
                                   title='', xlabel='', ylabel='',
                                   colors=['red', 'green', 'blue'],
                                   labels=['red', 'green', 'blue'],
                                   linestyle=['-', '-', '-', '-'],
                                   marker=['.', '.', '.', '.'],
                                   legend_fontsize=18, xlabel_fontsize=20,
                                   ylabel_fontsize=20,
                                   figsize=(6,5),
                                   legend_loc='upper right'):
    """Plot multiple results.
    
    NOTE: The time for each result should be along the time axis.
    """
    
    assert len(likelihood_arr) <= len(colors), "Missing colors for plot."
    mean_likelihood, std_likelihood = [], []

    if time_axis >= 0:
        for l in likelihood_arr:
            mean_likelihood.append(np.mean(l, axis=time_axis))
            std_likelihood.append(np.std(l, axis=time_axis))
    else:
        assert len(likelihood_arr[0].shape) == 1, \
            "Can only plot vector with -ve time axis"
        mean_likelihood = likelihood_arr
        std_likelihood = [np.zeros(likelihood_arr[0].shape)] * len(likelihood_arr)
        
    if x is None:
        x = range(len(mean_likelihood[0]))
    
    plots = []
    fig = plt.figure(figsize=figsize)
    ax = fig.add_subplot(111)
    
    for i in xrange(len(mean_likelihood)):
        m, s = mean_likelihood[i], std_likelihood[i]
        c = colors[i]
        p, = ax.plot(x, m, 'k', color=c, label=labels[i], 
                      linestyle=linestyle[i], 
                      # marker=marker[i],
                     )

        
        ax.fill_between(x, m - s, m + s,
                         alpha=0.2, edgecolor=c, facecolor=c,                         
                         linewidth=4, linestyle='dashdot', antialiased=True)
        

        plots.append(p)
        
    plt.legend(handles=plots, fontsize=legend_fontsize, loc=legend_loc)
    ax.set_title(title)
    ax.set_xlabel(xlabel, fontsize=xlabel_fontsize)
    ax.set_ylabel(ylabel, fontsize=ylabel_fontsize)

    if len(save_path) > 0:
        plt.savefig(save_path, bbox_inches="tight")
    plt.show()

    def plotPath(self, seq, poses_gt, poses_result):
        plot_keys = ["Ground Truth", "Ours"]
        fontsize_ = 20
        plot_num = -1

        poses_dict = {}
        poses_dict["Ground Truth"] = poses_gt
        poses_dict["Ours"] = poses_result

        fig = plt.figure()
        ax = plt.gca()
        ax.set_aspect('equal')

        for key in plot_keys:
            pos_xz = []
            # for pose in poses_dict[key]:
            for frame_idx in sorted(poses_dict[key].keys()):
                pose = poses_dict[key][frame_idx]
                pos_xz.append([pose[0, 3], pose[2, 3]])
            pos_xz = np.asarray(pos_xz)
            plt.plot(pos_xz[:, 0], pos_xz[:, 1], label=key)

        plt.legend(loc="upper right", prop={'size': fontsize_})
        plt.xticks(fontsize=fontsize_)
        plt.yticks(fontsize=fontsize_)
        plt.xlabel('x (m)', fontsize=fontsize_)
        plt.ylabel('z (m)', fontsize=fontsize_)
        fig.set_size_inches(10, 10)
        png_title = "sequence_{:02}".format(seq)
        plt.savefig(self.plot_path_dir + "/" + png_title + ".pdf", bbox_inches='tight', pad_inches=0)

def make_histograms(truescores, fpscores, title):
	plt.hist(truescores, normed=True, histtype='step', linewidth=3, label="True Events, %d" % len(truescores))	
	plt.hist(fpscores, normed=True, histtype='step', linewidth=3, label="FP Events, %d" % len(fpscores))
	plt.legend()
	plt.xlabel("Likelihood Score")
	plt.ylabel("Density of events")
	plt.title("Histogram of TP and FP Event Likelihoods")	

def plot_robots_ratio_time_micmac(deploy_robots_mic, deploy_robots_mac, deploy_robots_desired, delta_t):
    plot_option = 0 # 0: ratio, 1: cost
    num_iter = deploy_robots_mic.shape[1]
    total_num_robots = np.sum(deploy_robots_mic[:,0,:])
    
    diffmic_sqs = np.zeros(num_iter)
    diffmac_sqs = np.zeros(num_iter)
    diffmic_rat = np.zeros(num_iter)
    diffmac_rat = np.zeros(num_iter)
    for t in range(num_iter):
        diffmic = np.abs(deploy_robots_mic[:,t,:] - deploy_robots_desired)    
        diffmac = np.abs(deploy_robots_mac[:,t,:] - deploy_robots_desired) 
        diffmic_rat[t] = np.sum(diffmic) / total_num_robots       
        diffmic_sqs[t] = np.sum(np.square(diffmic))
        diffmac_rat[t] = np.sum(diffmac) / total_num_robots 
        diffmac_sqs[t] = np.sum(np.square(diffmac))
        
    x = np.arange(0, num_iter) * delta_t
    if(plot_option==0):
        l1 = plt.plot(x,diffmic_rat)
        l2 = plt.plot(x,diffmac_rat)
    if(plot_option==1):
        l1 = plt.plot(x,diffmic_sqs)
        l2 = plt.plot(x,diffmac_sqs)
    
    plt.xlabel('time [s]')    
    plt.ylabel('ratio of misplaced robots')
    plt.legend((l1, l2),('Micro','Macro'))
    plt.show()

 def _plot_histogram(self, data, number_of_devices=1, 
         preamp_timeout=1253):
     if number_of_devices == 0:
         return
     data = np.array(data)
     plt.figure(3)
     plt.ioff()
     plt.get_current_fig_manager().window.wm_geometry("800x550+700+25")
     plt.clf()
     if number_of_devices == 1: 
         plt.hist(data[0,:], bins=preamp_timeout, range=(1, preamp_timeout-1),
             color='b')
     elif number_of_devices == 2:
         plt.hist(data[0,:], bins=preamp_timeout, range=(1, preamp_timeout-1),
             color='r', label='JPM A')
         plt.hist(data[1,:], bins=preamp_timeout, range=(1, preamp_timeout-1),
             color='b', label='JPM B')
         plt.legend()
     elif number_of_devices > 2:
         raise Exception('Histogram plotting for more than two ' +
         'devices is not implemented.')
     plt.xlabel('Timing Information [Preamp Time Counts]')
     plt.ylabel('Counts')
     plt.xlim(0, preamp_timeout)
     plt.draw()
     plt.pause(0.05)

def visualize(u1, t1, u2, t2, U, omega):
    plt.figure(1)
    plt.plot(t1, u1, 'r--o')
    t_fine = np.linspace(0, t1[-1], 1001)  # мелкая сетка для точного решения
    u_e = u_exact(t_fine, U, omega)
    plt.hold('on')
    plt.plot(t_fine, u_e, 'b-')
    plt.legend([u'приближенное', u'точное'], loc='upper left')
    plt.xlabel('$t$')
    plt.ylabel('$u$')
    tau = t1[1] - t1[0]
    plt.title('$\\tau = $ %g' % tau)
    umin = 1.2*u1.min();
    umax = -umin
    plt.axis([t1[0], t1[-1], umin, umax])
    plt.savefig('tmp1.png');  plt.savefig('tmp1.pdf')
    plt.figure(2)
    plt.plot(t2, u2, 'r--o')
    t_fine = np.linspace(0, t2[-1], 1001)  # мелкая сетка для точного решения
    u_e = u_exact(t_fine, U, omega)
    plt.hold('on')
    plt.plot(t_fine, u_e, 'b-')
    plt.legend([u'приближенное', u'точное'], loc='upper left')
    plt.xlabel('$t$')
    plt.ylabel('$u$')
    tau = t2[1] - t2[0]
    plt.title('$\\tau = $ %g' % tau)
    umin = 1.2 * u2.min();
    umax = -umin
    plt.axis([t2[0], t2[-1], umin, umax])
    plt.savefig('tmp2.png');
    plt.savefig('tmp2.pdf')

def visualize(segmentation, expression, visualize=None, store=None, title=None, legend=False):
    notes = []
    onsets = []
    values = []
    param = ['Dynamics', 'Articulation', 'Tempo']
    converter = NoteList()
    converter.bpm = 100
    if not visualize:
        visualize = selectSubset(param)
    for segment, expr in zip(segmentation, expression):
        for note in segment:
            onsets.append(converter.ticks_to_milliseconds(note.on)/1000.0)
            values.append([expr[i] for i in visualize])
    import matplotlib.pyplot as plt
    fig = plt.figure(figsize=(12, 4))
    for i in visualize:
        plt.plot(onsets, [v[i] for v in values], label=param[i])
    plt.ylabel('Deviation')
    plt.xlabel('Score time (seconds)')
    if legend:
        plt.legend(bbox_to_anchor=(0., 1), loc=2, borderaxespad=0.)

    if title:
        plt.title(title)
    #dplot = fig.add_subplot(111)
    #sodplot = fig.add_subplot(111)
    #dplot.plot([i for i in range(len(deltas[0]))], deltas[0])
    #sodplot.plot([i for i in range(len(sodeltas[0]))], sodeltas[0])
    if store:
        fig.savefig('plots/{0}.png'.format(store))
    else:
        plt.show()

    def make_bar(
        x,
        y,
        f_name,
        title=None,
        legend=None,
        x_label=None,
        y_label=None,
        x_ticks=None,
        y_ticks=None,
    ):
        fig = plt.figure()

        if title is not None:
            plt.title(title, fontsize=16)
        if x_label is not None:
            plt.ylabel(x_label)
        if y_label is not None:
            plt.xlabel(y_label)
        if x_ticks is not None:
            plt.xticks(x, x_ticks)
        if y_ticks is not None:
            plt.yticks(y_ticks)

        plt.bar(x, y, align="center")

        if legend is not None:
            plt.legend(legend)

        plt.savefig(f_name)
        plt.close(fig)

def plotISVar():
    plt.figure()
    plt.title('Variance minimization problem (call).\nVertical lines mark the minima.')
    for K in [0.6, 0.8, 1.0, 1.2]:
        theta = np.linspace(-0.6, 2)
        var = [BS.exactCallVar(K*s0, theta) for theta in theta]
        minth = theta[np.argmin(var)]
        line, = plt.plot(theta, var, label=str(K))
        plt.axvline(minth, color=line.get_color())

    plt.xlabel(r'$\theta$')
    plt.ylabel('call variance')
    plt.legend(title=r'$K/s_0$', loc='upper left')
    plt.autoscale(tight=True)

    plt.figure()
    plt.title('Variance minimization problem (put).\nVertical lines mark the minima.')
    for K in [0.8, 1.0, 1.2, 1.4]:
        theta = np.linspace(-2, 0.5)
        var = [BS.exactPutVar(K*s0, theta) for theta in theta]
        minth = theta[np.argmin(var)]
        line, = plt.plot(theta, var, label=str(K))
        plt.axvline(minth, color=line.get_color())

    plt.xlabel(r'$\theta$')
    plt.ylabel('put variance')
    plt.legend(title=r'$K/s_0$', loc='upper left')
    plt.autoscale(tight=True)

    def make_line(
        x,
        y,
        f_name,
        title=None,
        legend=None,
        x_label=None,
        y_label=None,
        x_ticks=None,
        y_ticks=None,
    ):
        fig = plt.figure()

        if title is not None:
            plt.title(title, fontsize=16)
        if x_label is not None:
            plt.ylabel(x_label)
        if y_label is not None:
            plt.xlabel(y_label)
        if x_ticks is not None:
            plt.xticks(x, x_ticks)
        if y_ticks is not None:
            plt.yticks(y_ticks)

        if isinstance(y[0], list):
            for data in y:
                plt.plot(x, data)
        else:
            plt.plot(x, y)

        if legend is not None:
            plt.legend(legend)

        plt.savefig(f_name)
        plt.close(fig)